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FOREWORD 


This  computer  program  User's  Manual  was  prepared 
by  the  Los  Angeles  Division  of  Science  Applications,  In- 
corporated, Los  Angeles,  California  for  the  Vehicle  Dy- 
namics Division  of  the  Air  Force  Flight  Dynamics  Labora- 
tory, Wright-Patterson  Air  Force  Base,  Ohio.  The  computer 
programs  were  developed  under  Project  1370,  "Dynamic  Prob- 
lems in  Flight  Vehicles,"  Task  137004,  "Design  Analysis," 
Contract  F33615-74-C-3094 . James  J.  Olsen  and  later  Lt. 
William  L.  Holman  of  AFFDL/FYS  were  the  Air  Force  task 
engineers . 

R.  M.  Traci  was  the  principal  investigator  for 
the  study  and  J.  L.  Farr,  Jr.  developed  the  computer  pro- 
grams described  in  this  report.  Consultant  E.  D.  Albano 
contributed  to  the  development  and  implementation  of  the 
numerical  method. 

The  authors  submitted  this  report  in  October  1974 
for  publication  as  an  AFFDL  technical  report  to  cover  re- 
search performed  from  June  through  September  1974. 


-iii- 


TABLE  OF  CONTENTS 


Page 

1.0  INTRODUCTION  1 

2 . 0 SMALL  PERTURBATION  THEORY  FOR  UNSTEADY  3 

TRANSONIC  FLOW 

3.0  SUMMARY  OF  NUMERICAL  SOLUTION  PROCEDURE  8 

3.1  Numerical  Solution  Procedure  8 

for  the  Steady  Potential 

3.1.1  Finite  Difference  Forms  9 

and  Iterative  Solution 
Procedure 

3.1.2  Airfoil  Boundary  Condition  13 

3.1.3  Kutta  Condition  15 

3.1.4  Farfield  Boundary  Condition  18 

3.2  Numerical  Solution  Method  for  19 

First  Order  Unsteady  Perturbation 
Potential 

4.0  PROGRAM  DESCRIPTIONS  25 

5.0  INPUT  AND  OUTPUT  32 

5.1  STRANS  Input  34 

5.2  STRANS  Output  38 

5.3  UTRANS  Input  39 

5.4  UTRANS  Output  43 

-i  V- 


TABLE  OF  CONTENTS  (Cont'd) 


Page 

6.0  PROGRAM  USAGE  45 

6.1  Computational  Details  45 

6.2  Summary  of  Recommended  Compu-  50 

tational  Procedures 

7.0  SAMPLE  CASES  52 

7.1  STRANS  Test  Cases  52 

7.1.1  Input  for  STRANS  Test  Cases  54 

7.1.2  Sample  Output  for  STRANS  59 

Test  Cases 

7.2  UTRANS  Test  Cases  71 

7.2.1  Input  for  UTRANS  Test  Cases  73 

7.2.2  Sample  Output  for  UTRANS  78 

Test  Cases 

8.0  REFERENCES  90 

APPENDIX  A:  FORTRAN  LISTING  OF  STRANS  91 

APPENDIX  B:  FORTRAN  LISTING  OF  UTRANS  108 


-V- 


LIST  OF  ILLUSTRATIONS 


FIGURE  1 

FIGURE  2 

FIGURE  3 

FIGURE  4 
FIGURE  5 

FIGURE  6 
FIGURE  7 


SCHEMATIC  OF  AIRFOIL  GEOMETRY 
AND  TRANSONIC  FLOWFIELD 


SCHEMATIC  OF  NUMERICAL  SOLUTION 
DOMAIN 


SKETCH  OF  LOCAL  MESH  ABOUT  A 
COMPUTATIONAL  POINT 


SCHEMATIC  OF  MESH  NEAR  AIRFOIL 


SCHEMATIC  OF  NUMERICAL  SOLUTION 
DOMAIN 


LOGICAL  FLOW  OF  STRANS  PROGRAM 


DEFINITION  OF  PRESSURE  AND  FORCE 
COEFFICIENTS 


Page 


3 


9 


10 


13 

20 


26 

33 


-vi- 


< M 


LIST  OF  SYMBOLS 


a 

Sound  speed,  coefficients  in  asymp- 
totic expansion 

c 

Airfoil  chord 

C , C 
P P 

Unsealed  and  scaled  pressure  coeffi- 
cient 

ACP 

Jump  in  pressure  coefficient  over 
airfoil 

Cl' 

Unsealed  and  scaled  lift  coefficient 

EPSC0L, EPSGRD 

Variables  to  control  column  and  grid 
convergence 

f,t 

Airfoil  shape  function,  thickness  dis- 
tribution 

«(2) 

Hankel  function  of  second  kind 

i 

Designates  imaginary  part 

i»  j 

Mesh  indices 

I , I 

1 2 

Integrals  in  farfield  expression 

k 

A 

k 

Reduced  frequency  based  on  chord 

Ratio  of  scaled  frequency  to  transonic 
similarity  parameter 

K 

Transonic  similarity  paramter 

M 

OO 

Freestream  Mach  number 

M, NGFF , PGFF 

Parameters  in  grid  and  circulation 
iteration  schemes 

R,  Rx,  S 

Transformed  coordinates 

<\j 

t , t 

Physical  and  scaled  time 

U 

Freestream  velocity 

-vi  i - 


LIST  OF  SYMBOLS  (Cont'd) 


V 

Coefficient  in  transonic  potential 
equation 

% 

X,  X 

Physical  and  scaled  streamwise  coordi- 
nate 

'Xj 

Y'Y 

Physical  and  scaled  normal  coordinate 

a 

Angle  of  attack 

6 

Airfoil  thickness  parameter 

A 

Mesh  spacing  or  jump  in  quantity 

e 

Unsteady  perturbation  parameter 

Y 

Ratio  of  specific  heats 

Y , Y 
0 1 

Steady  and  unsteady  circulation 

U)  , 

Physical  and  scaled  frequency  of  os- 
cillation 

0)  . 
D 

Relaxation  parameter 

a 

Jump  in  potential  across  wake 

<P  r$ 

Unsealed  and  scaled  perturbation 
velocity  potential 

V 

Total  velocity  potential 

Subscripts 


oo 

Refers  to  freestream 

if  j 

Refers  to  finite  difference  grid  points 

ff 

Refers  to  farfield 

te 

Refers  to  trailing  edge  of  airfoil 

a,  l 

Refers  to  upper  and  lower  airfoil  sur- 
face 

x,  y , t 

Refers  to  partial  derivatives  with 
respect  to  coordinates  and  time 

LIST  OF  SYMBOLS  (Cont'd) 


Superscripts 


0 r 1 

Refers  to  zero  order  steady  and  first 
order  unsteady  perturbations 

E 

Refers  to  elliptic,  centered  differ- 
ence form 

H 

Refers  to  hyperbolic,  backward  differ 
ence  form 

V 

Refers  to  iteration  number 

+ / - 

Refers  to  above  or  below  airfoil  or 
wake 

-ix- 


1.0 


INTRODUCTION 


Computer  programs  STRANS  and  UTRANS  implement  a 
small  disturbance  potential  flow  theory  for  the  two- 
dimensional  unsteady  transonic  flow  about  thin  airfoils 
undergoing  low  reduced  frequency  harmonic  oscillations. 

The  theory  is  based  on  Landahl ' s suggestion1'2  that  a 
linear  system  can  be  obtained  by  considering  the  unsteady 
flow  as  a small  perturbation  to  the  mean  nonuniform  flow. 
The  perturbation  expansion  approach  has  recently  been  de- 
veloped with  different  emphasis  in  independent  studies  by 
Traci,  et  al.3,  and  Ehlers4 . The  computer  programs,  docu- 
mented in  this  users  manual,  were  developed  during  the 
former  study  and  detailed  descriptions  of  the  theory  and 
numerical  solution  technique  with  calculated  results  are 
presented  in  the  final  report  of  that  study3. 

In  the  perturbation  expansion  approach  used,  the 
perturbation  potential  function  is  expanded  in  a series  of 
increasing  powers  of  a small  parameter  which  is  a measure 
of  the  amplitude  of  an  unsteady  disturbance  to  the  bound- 
ary. The  resulting  expansion  of  the  unsteady  potential 
equation  results  in  a sequence  of  partial  differential 
equations  for  the  perturbation  potentials.  The  zeroth 
order  equation  is  the  usual  nonlinear  steady  transonic  po- 
tential equation  of  mixed  elliptic/hyperbolic  type  and  is 
solved  in  STRANS  using  the  mixed  differencing,  relaxation 
procedure  of  Murman  and  Cole4.  The  first  order  unsteady 
potential  equation  is  linear  and  for  harmonic  boundary  dis- 
turbances is  also  of  the  mixed  elliptic/hyperbolic  type, 
depending  upon  the  steady  solution.  It  is  solved  in  UTRANS 
using  the  same  numerical  technique  as  used  in  STRANS. 

The  background,  approximations  and  practical  accu- 
racy of  the  theory  and  numerical  solution  procedure  are 
discussed  in  some  detail  in  the  companion  report3  to  this 
users  manual.  As  noted  in  that  report,  the  method  imple- 
mented in  the  present  versions  of  the  programs  has  not 
been  completely  verified  as  to  its  accuracy  due  to  the 
dearth  of  reliable  experimental  data  or  accurate  alternate 
solutions  for  comparison.  This  lack  of  empirical  or  ana- 
lytical techniques  for  predicting  unsteady  aerodynamic 
forces  in  the  transonic  speed  range  is  in  fact  the  justi- 
fication for  publishing  the  preliminary  versions  of  the 
computer  programs  described  here.  The  approach  embodied 
in  the  programs  shows  considerable  promise  of  providing  a 
practical  method  for  filling  this  void.  There  are  many 
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ways  theoretically,  numerically  and  operationally  that 
the  present  versions  of  STRANS  and  UTRANS  may  be  improved. 
Work  is  continuing  in  this  direction,  so  that  it  is  antic- 
ipated that  the  efficiency  and  capability  of  the  programs 
will  continue  to  increase  in  the  future. 

The  theory  and  practice  of  the  computer  program 
operation  are  discussed  in  the  following  sections.  The 
small  perturbation  theory  and  numerical  solution  proce- 
dure are  summarized  in  Sections  2.0  and  3.0,  respectively. 

A description  of  the  program's  logical  operation  and  a 
brief  subroutine  description  are  given  in  Section  4.0. 

Section  5.0  presents  a complete  description  of  the  program 
input,  with  suggested  values  for  various  control  variables, 
and  the  program  output.  Section  6.0  describes  the  program 
usage  and  includes  suggestions  for  making  effective  use  of 
the  programs.  Past  experience  with  convergency,  accuracy 
and  efficiency  is  discussed  as  well  as  a summary  of  a recom- 
mended computational  procedure.  Sample  cases  which  exercise 
all  program  options  are  presented  in  Section  7.0  with  a com- 
plete specification  of  all  input  and  sample  output.  Finally, 
complete  FORTRAN  listings  of  STRANS  and  UTRANS  are  pre- 
sented in  the  appendices. 
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2.0 


SMALL  PERTURBATION  THEORY  FOR  UNSTEADY  TRANSONIC 
FLOW 


Small  disturbance  theory  is  the  principal  analyti- 
cal tool  for  all  speed  ranges  and  the  transonic  speed 
range  is  no  exception.  The  transonic  flow  regime  (de- 
fined roughly  as  | 1 - M^ | ^6  where  M is  the  freestream 
Mach  number  and  6 is  the  airfoil  thickness  ratio)  is  in- 
herently nonlinear  however,  so  that  linear  theories  which 
apply  to  fully  subsonic  or  supersonic  flow  are  not  valid. 
Special  account  must  therefore  be  taken  of  the  nonlinear- 
ity which,  of  course,  severely  complicates  the  formulation 
so  that  finite  difference  solution  procedures  are  indicated. 
The  theory  summarized  here  is  based  upon  the  treatment  of 
the  unsteady  flow  as  a small  perturbation  of  the  steady 
transonic  flow.  This  effectively  linearizes  the  oscillat- 
ing component  of  the  flow  about  the  nonlinear  steady  solu- 
tion. The  resulting  boundary  value  problems  for  the  steady 
(solved  in  STRANS)  and  first  order  unsteady  (solved  in 
UTRANS)  perturbation  potentials  are  presented  in  this  sec- 
tion. 


Figure  1.  Schematic  of  Airfoil  Geometry  and 
Transonic  Flowfield 
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The  problem  of  interest  is  the  two-dimensional 
flow  about  an  airfoil  oscillating  in  shape,  angle  of  at- 
tack or  flap  angle  in  the  transonic  flow  speed  range.  Of 
particular  interest  are  the  steady  and  unsteady  aerodynam- 
ic forces  on  the  airfoil  for  use  in  flutter  studies.  The 
airfoil  geometry,  flowfield  schematic  and  coordinate  defi- 
nition are  given  in  Figure  1 above.  Rectangular  coordi- 
nates (x,y)  are  fixed  to  the  airfoil  leading  edge  and  U, 

M^,  a^  are  the  freestream  velocity,  Mach  number  and  sound 
speed  respectively.  The  airfoil  has  a thickness  ratio  6 
which  is  the  airfoil  maximum  thickness  divided  by  its  chord  c 
and  an  angle  of  attack  a.  The  assumption  is  made  that 
6 <<  1 and  a is  of  the  same  order  of  magnitude  as  6.  Al- 
so the  oscillatory  motion  of  the  airfoil  is  assumed  to  be 
described  by  a small  non-dimensional  displacement  e <<  6 
and  a reduced  frequency  k = wc/U  based  on  airfoil  chord 
where  w is  the  frequency  of  oscillation.  The  reduced  fre- 
quency range  treated  by  the  present  analysis  in  k <<  1 or 
more  explicitly  k of  order  62/3. 

Assuming  inviscid,  isentropic  flow,  the  problem  can 
be  reduced  to  the  solution  of  a single  equation  for  a ve- 
locity potential  plus  the  tangency  boundary  condition  on 
the  airfoil  surface.  As  is  well  known  the  derivation  of 
a small  disturbance  theory  for  transonic  flows  requires  a 
singular  perturbation  approach.  For  the  low  frequency  un- 
steady flow  regime  of  interest,  the  following  scaling  is 
introduced : 
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and  the  total  potential  is  expanded  about  the  uniform  flow 
thusly : 
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Retaining  all  terms  to  leading  order  in  the  total  poten- 
tial equation  and  boundary  condition  results  in  the  fol- 
lowing form  for  the  unsteady  small  perturbation  system: 
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where  the  transonic  similarity  parameter  is: 
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where  f„  o is  the  unsteady  airfoil  shape  function  (Equa- 
tion 7 below)  on  the  upper  and  lower  surface^  respectively, 
and  where  []  denotes  a jump  between  0 and  0 . It  is 
noted  that  the  airfoil  tangency  boundary  condition  (Equa- 
tion 4)  and  Kutta  condition  (Equation  5)  are  applied  in 
the  small  disturbance  manner  on  y = 0 and  in  the  low  fre- 
quency, quasi-steady  sense  (k  -*•  0)  . 

The  system  of  Equations  3 through  6 provide  a for- 
mulation of  the  unsteady  airfoil  problem  in  the  nonlinear 
domain,  which  includes  flowfields  with  shocks.  The  formu- 
lation is  self  consistent  for  similarity  parameter  K of 
order  1 and  for  reduced  frequency  k of  order  62/3. 

The  approach  to  solving  the  nonlinear  system  given 
above  (Equations  3 through  6)  is  to  expand  the  perturba- 
tion potential  function  in  terms  of  the  unsteady  boundary 
disturbance  e <<  1.  From  this  point  on  all  tildas  (M 
will  be  dropped  with  the  understanding  that  all  variables 
are  scaled  variables.  Harmonic  boundary  disturbances  are 
explicitly  treated: 

'u.l1*'*  = fi,£(x)  + e fi,£<x>eint  <7> 

where  the  scaled  frequency  is: 
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The  perturbation  potential  is  expanded  as  follows: 

<Mx,y,t)  = <J>°(x,y)  + e <j> 1 (x,  y ) eifit  + . . . (9) 


Substituting  this  into  the  perturbation  potential  equation 
plus  boundary  conditions  and  combining  terms  results  in 
the  following  pair  of  boundary  value  problems  for  <p 0 and 
<t> 1 respectively. 


and 


(K-(j) 0 ) (j)  0 + <p°  = 0 

x Txx  Tyy 


<py  = f“;yx),  on  y = ±0,  0 ;<  x £ 1 


[<J>°]  = 0,  on  y = 0,  x > 1 

X 


( ” ) 2 + ( 4>  ° ) 2 -»■  o,  as  x2  + y2 
a y 


(K^x^ix  + ^yy  - (Cc  + 2ifi)^  = 0 


(10) 


II 

-e- 

fu',l 

(x)  , 

ony=±0,  0<x 

= 0, 

on  y = 0,  x > 1 

-e- 

2 + 

( ^'y ) 2 

0,  as  x2  + y2 

(ID 


Equation  10  is  recognized  as  the  usual  formulation  for 
steady  transonic  flow  and  Equation  11  is  the  formulation 
for  the  low  frequency  unsteady  perturbation  thereof.  Note 
that  the  governing  equation  for  <f) 1 is  linear  but  of  the 
same  mixed  elliptic/hyperbolic  type  as  the  steady  poten- 
tial equation  depending  upon  the  steady  solution.  It  is 
also  noted  that  cf> 1 is  in  general  complex  thereby  permitting 
phase  shifts  between  the  field  quantities  and  the  boundary 
disturbance.  A more  general  formulation  including  higher 
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order  unsteady  modes  and  description  of  the  nature  of  the 
equations  is  given  in  Reference  3. 

The  main  physical  quantities  of  interest  are  the 
pressure  coefficient  and  airfoil  force  coefficients.  The 
pressure  coefficient,  defined  in  the  usual  manner,  is  giv- 
en by: 


6 2/ 3 

[ (l+y)M2] l/3 


(C°  + 
P 


iftt. 
e ) 


(12) 


where  the  steady  and  unsteady  scaled  pressure  coefficients 
are  given  to  leading  order  in  the  small  disturbance  and 
low  frequency  approximation  by: 


Also  the  lift  coefficient  is  given  by: 


(13) 


C 


Z 


26  2/ 3 

[ ( 1+y ) M^] !/3 


(Y 

o 


+ e 


Y e 
i 


iftt  j 


(14) 


where  y„  = [<J>  1 i and  yx  = 1 1 ] x:=  l are  the  steady  and  un- 
steady circulations  about  the_airfoil.  The  finite  differ- 
ence forms  used  to  calculate  Cp  and  Cp  are  defined  in  the 
next  section  and  the  detailed  definition  of  the  calculated 
generalized  force  coefficients  is  presented  in  Section  4.0. 
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3.0 


SUMMARY  OF  NUMERICAL  SOLUTION  PROCEDURE 


The  numerical  solution  procedures  for  the  boundary 
value  problems  for  the  steady  perturbation  potential  (Equa- 
tion 10)  and  the  unsteady  perturbation  potential  (Equation 
11) , as  implemented  in  computer  programs  STRANS  and  UTRANS 
respectively,  are  now  summarized.  As  has  been  pointed  out, 
both  equations  are  of  the  same  mixed  elliptic/hyperbolic 
type  and  are  solved  using  the  mixed  differencing  relaxa- 
tion procedure  of  Murman  and  Cole5.  Summary  of  the  pro- 
cedure used  in  STRANS  for  calculating  the  steady  solution 
is  given  in  Section  3.1.  Variations  on  the  basic  procedure, 
used  in  UTRANS  for  the  unsteady  perturbation  are  discussed 
in  Section  3.2.  It  is  noted  that  the  finite  difference 
equations  which  follow  are  written  in  terms  of  the  actual 
FORTRAN  variables  used  in  the  programs  to  as  great  a de- 
gree as  possible.  In  addition,  following  usual  program- 
ming practice,  the  FORTRAN  variables  are  descriptive  of 
the  physical  variables. 


3 . 1 Numerical  Solution  Procedure  for  the  Steady  Potential 


The  difference  scheme  and  relaxation  solution  pro- 
cedure summarized  in  this  section  are  based  on  the  work 
of  Murman5'6'8,  Cole5  and  Krupp6'7.  They  pointed  out  the 
essential  ingredient  for  the  success  of  relaxation  proce- 
dures for  the  steady  transonic  potential  equation.  The  key 
to  their  approach  is  to  account  for  the  local  nature  of 
the  flow  (elliptic  in  subsonic  regions,  hyperbolic  in  su- 
personic regions)  in  the  finite  difference  approximation 
to  the  governing  equations.  The  version  of  this  technique, 
as  implemented  in  STRANS,  is  summarized  in  this  section.  A 
more  detailed  description  can  be  found  in  Reference  3. 

Consider  the  steady,  transonic  potential  equation, 
with  boundary  conditions,  in  the  numerical  solution  domain 
indicated  schematically  below.  Finite  difference  grid  lines 
are  defined  in  the  region  around  the  airfoil  which  is  fixed 
on  the  slit  0 <_  x <_  1,  y = 0.  Values  of  the  potential  are 
defined  at  the  intersection  of  the  grid  lines  (i,j).  The 
finite  difference  approximations  to  the  governing  equation 
at  a general  grid  point  account  for  the  varying  type  of 
the  equation  as  is  summarized  in  Section  3.1.1.  An  itera- 
tive line  relaxation  solution  procedure,  also  described  in 
Section  3.1.1,  uses  the  finite  difference  equations  to 
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successively  update  each  point  in  the  grid  starting  from 
some  initial  guess.  The  special  numerical  treatments  re- 
quired for  the  airfoil  boundary,  Kutta  condition  and  far- 
field  are  described  in  Sections  3.1.2  to  3.1.4,  respec- 
tively. 


j = JM 


Figure  2,  Schematic  of  Numerical  Solution  Domain 


3.1.1  Finite  Difference  Forms  and  Iterative  Solution 
Procedure 


In  the  numerical  scheme,  a rectangular  mesh  with 
general  grid  line  spacing,  as  indicated  in  the  sketch  (Fig- 
ure 3) , is  used.  Uneven  grid  spacing  makes  it  possible 
to  concentrate  grid  points  near  the  airfoil  and  in  regions 
where  rapid  changes  in  the  potential  or  its  derivative 
(wing  leading  edge,  shocks,  etc.)  are  expected.  Expansion 
of  the  grid  spacing  away  from  the  airfoil  out  to  the 
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boundaries  of  the  mesh  (i  = 1 or  IM,  j = 1 or  JM)  per- 
mits economic  use  of  grid  points  while  maximizing  resolu- 
tion in  the  region  of  interest. 


Y>  j 


x,  i 


Ay: 

D , 

Ax . ~ 
1-2 

Ax . . 

l- 1 

Ax  . 

l 

Ay:-i 

Figure  3.  Sketch  of  Local  Mesh  About  A 
Computational  Point 


At  a general  point  (i,j)  in  the  finite  difference 
mesh,  the  solution  begins  by  center  differencing  the  co- 
efficient of  in  the  following  manner: 


V. 

1 


j 


(K-4>°).  . 

Tx'i,j 


= K 


- AX1  . 

l 


( <f>  . . , . -<p  . .) 

yi+l/ 1 i,D 


- AX2  i 


(15) 


with 


Ax . . 

AX1 . = ^ 1 

i Ax  . (Ax  . 


. . +Ax . ) 
l-l  l 


AX2  . 
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Ax. 

l 

Ax . (Ax . . +Ax . ) 
l-l  l-l  l 


-10- 


Tests  are  then  made  on  this  coefficient  to  determine  the 
nature  of  the  potential  equation  at  the  point  (i,j).  if 
vi  j > 0 the  equation  is  elliptic  so  that  <J>£X  is  center 
differenced  and  if  Vj^j  < 0 the  equation  is  hyperbolic 
and  backward  differencing  is  used  for  both  Vi  t j and  (j>£x. 

An  additional  test  is  made  to  determine  if  the  point  is 
near  the  sonic  line.  If  Vi,j  < 0 but  > 0,  the 

point  is  identified  as  a "parabolic  poinf’and  the  center 
differenced  form  of  V±r j is  used  but  backward  differenced 
form  for  <j>xx  is  used.  The  4>yy  term  is  always  center  dif- 
differenced  regardless  of  the  above  considerations.  The 
finite  difference  form  of  the  transonic  potential  equation 
can  thus  be  summarized  as  follows: 


V 


• . > 0,  elliptic  (subsonic) 

±-LJ ■ 


V.  . (4)  )?  . + (d>  ) . . = 0 

i,]  xx'i,]  yy  1,1 

(16) 

V.  . < 0, 

-J-J 

hyperbolic  (supersonic) 

v1?  . ( <})  )h.  + (4>  ).  . = 0 

i,]  xx'i,]  yy  i,i 

(17) 

V.  . < 0, 
-JLt J 1 

V.  . . > 0,  parabolic  (sonic) 

i i ' 1 

Vi,j(l{>xx)i,j  “ ^yy^  i,  j ~ 0 

(18) 

where : 

VH  . = K 
i/l 

(19) 

<*xx'i 

(20) 

(21) 
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(cp  ) . . = AYl  . (4>  . . .. -<p  . .)  - AY 2 . (4)  . .-4)  . . . ) (22) 

yy  i/]  3 i/3+i  *1,3'  *1,3-1' 


and  where 


BXl  . = 2 AXl./Ax.  . , BX2.  = 2 AX2 . /Ax . 
1 i'  i-l'  i r l 


CX.  = 1/ ( 2Ax . ) 

l i 


(23) 


AYlj  AyTTAyT+AyT^^  ' AY2j  Ayj_1  (Ay^+Ay..^) 


Using  the  above  forms,  the  finite  difference  equa- 
tions are  set  up  for  each  column  (x  = constant)  in  turn, 
taking  into  account  airfoil,  wake  and  farfield  boundary 
conditions.  This  results  in  a sequence  of  nonlinear  alge- 
braic equations  for  each  column  of  4>'s,  which  can  be  solved 
by  iteration  or  Newton's  method.  The  method  used  in 
STRANS  is  to  linearize  the  system  of  equations  by  using 
the  previous  iterate  for  the  coefficient  Vj^j  = K - 4>x- 
The  resulting  linear  system  is  tridiagonal  and  is  easily 
solved  by  Gaussian  elimination.  The  iteration  process  is 
terminated  when  the  difference  between  successive  iterates, 
for  all  points  in  the  column,  is  less  than  some  specified 
small  amount  (EPSC0L) . The  process  works  quite  well,  re- 
quiring but  three  or  four  iterations  to  converge  to  good 
accuracy  (A4>  < 10-5). 

After  each  column  is  solved,  it  is  relaxed  with  a 
variable  relaxation  factor  which  depends  on  the  local  na- 
ture of  the  solution. 


<tv+1  = 4,u  . 
i.J  i.] 


,-vV+l  ,v  > 

w . ( <p . . ~4>  • • ) 
3i/3  i/3 


(24) 


j = 1,2,  ...  JM 


where  w is  the  relaxation  factor  and  $ is  the  solution  of 
of  the  system  of  equations  for  a column.  Relaxation  fac- 
tors which  have  worked  well  are  to  = 1.5  ■*  1.7  for  elliptic 
points  and  to  = .75  for  hyperbolic  points. 

The  column  solution  process  described  above  is  per- 
formed for  each  column  in  turn,  starting  from  an  initial 
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guess.  The  initial  guess  typically  is  a previous  calcula- 
tion for  the  same  or  a similar  airfoil  at  conditions  of 
Mach  number  and  angle  of  attack  close  to  the  problem  of 
interest.  If  such  a solution  doesn't  exist,  the  farfield 
expression  is  used  as  an  initial  guess  with  an  arbitrary 
truncation  of  the  singularity  at  x = y = 0.  The  grid  is 
swept  in  the  above  manner  until  the  change  in  <J>  for  all 
grid  points  during  one  grid  sweep  is  less  than  some  arbi- 
trary small  amount  ( EPSGRD ) . Typically  convergence  is 
quite  rapid  initially  but  becomes  slower  as  the  iteration 
process  proceeds. 

A technique  which  considerably  aids  convergence  is 
to  begin  the  calculation  on  a relatively  coarse  grid  and 
to  periodically  refine  the  grid  as  the  solution  proceeds. 

The  rationale  behind  the  technique  is  that  the  coarse  grid 
speeds  transmission  of  information  (circulation,  boundary 
conditions,  etc.)  throughout  the  grid.  Successive  refine- 
ment of  the  grid  then  increases  solution  accuracy  with  the 
interpolated  values  of  <p  from  the  coarse  grid  providing  a 
good  starting  point.  Currently  the  technique  is  implemented 
by  a simple  halving  of  the  grid  spacing. 


3.1.2  Airfoil  Boundary  Condition 


Figure  4,  Schematic  of  Mesh  Near  Airfoil 
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As  derived  above,  the  airfoil  tangency  boundary  con- 
dition is  applied  on  the  line  y=0,  0 £ x < 1,  in  the 
usual  small  perturbation  sense,  rather  than  on  the  physical 
boundary  of  the  airfoil.  This  is  a significant  simplifi- 
cation relative  to  the  finite  difference  approximation 
since  it  precludes  the  possibility  of  accounting  for  the 
details  of  the  boundary  geometry  in  the  finite  difference 
mesh.  Thus  the  rectangular  grid  used  throughout  the  solu- 
tion field  is  also  appropriate  near  the  airfoil.  The  air- 
foil boundary  conditions  fit  naturally  into  the  finite  dif- 
ference equations  for  a column  of  grid  points  as  required 
by  the  line  relaxation  procedure.  The  way  that  this  is 
accomplished  is  now  described. 

The  schematic  given  above  defines  the  mesh  struc- 
ture near  the  slit  on  which  the  airfoil  boundary  condi- 
tions are  applied;  ILE  and  ITE  define  the  columns  at  the 
leading  edge  and  trailing  edge  respectively  and  JW  defines 
the  row  location  of  the  slit.  For  any  column  such  that 
ILE  £ i _<  ITE  the  airfoil  boundary  condition  effectively 
closes  the  system  of  equations  for  <J>  at  grid  points  above 
and  below  the  slit.  The  mathematical  statement  of  the 
boundary  condition  is  recalled: 


<t>  = f »(x)  ony=± 

Yy  a, l * 


= ±0,  0 < x < 1 


(25) 


where  f^'i  are  defined  by  the  user  in  the  functions  FUP  and 
FLP  respectively.  The  essence  of  the  method  is  to  use 
Equation  25  in  the  finite  difference  approximation  to 
the  c|>yy  term  of  the  potential  equation  at  points  just 
above  (j  = JW+1)  and  just  below  (j  = JW-1)  the  airfoil. 

The  form  taken  for  this  term  is  as  follows: 


(4>  ) • • 

YY  irJ 


{AyT+TFy 


<t>lf  - f 0 * (x.  )j  (26) 

1 j-l)  ( Ayj  U 1 ( 


j = JW+1,  ILE  < i < ITE 


(<f>  ) • • 

yy  i»] 


tf°»'(xj  - 


( 0 • .-<(>.  • , ) 

i,]  1,1-1 


4yj-i 


(iy^^Ay^  ''■l 

j = JW-1,  ILE  < i < ITE 


(27) 
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where  subscript  u,£  refer  to  upper  or  lower  airfoil  sur- 
face respectively.  These  finite  difference  forms  effec- 
tively fix  the  boundary  condition  on  the  slit;  y = 0+  for 
the  upper  surface  and  y = 0“  for  the  lower  surface.  They 
also  are  consistent  with  the  second  order  accuracy  of  the 
differencing  at  general  points  in  the  mesh.  The  other 
terms  in  the  potential  equation  are  treated  in  the  same 
manner  as  at  general  field  points. 

Values  of  potential  on  the  airfoil  surface  are  not 
explicitly  calculated  in  the  scheme  used  here.  These  val- 
ues are  obtained  only  after  the  equations  for  a column  of 
grid  points  have  been  solved  and  the  values  of  <f>l  relaxed. 
<j>  on  the  airfoil  surface  is  then  obtained  by  linear  extrap- 
olation from  above,  for  the  upper  surface,  and  from  below 
for  the  lower  surface.  The  pressure  coefficient  on  the  air- 
foils is  calculated,  using  the  extrapolated  values  of  po- 
tential, by  center  differencing  using  the  following  equa- 
tion: 


2AX1.  (<j>.  . .-d>  . 

1 yi+1,d 


) - 2 AX  2 . ((}). 


. — 4> . . 

/ 3 i-l, 


) (28) 


The  airfoil  circulation  Yfce  is  another  quantity  that  is 
of  interest  and  which  is  needed  for  the  treatment  of  the 
Kutta  condition  discussed  in  the  next  section.  This  is 
calculated  in  a straightforward  manner,  using  the  values 
of  potential  at  the  trailing  edge  as  follows: 


(Vte  = '♦"’x-l  = *u|  - *1|  , <29> 

1 X=1  1 X=1 

where,  as  before,  subscripts  u and  £ refer  to  the  extrap- 
olated values  on  the  upper  and  lower  airfoil  surfaces  re- 
spectively. 

The  generalized  forces  on  the  airfoil  are  calculated  in 
STRANS  using  a trapezoidal  rule  integration  over  the  airfoil 
of  the  difference  in  upper  and  lower  surface  pressure  coef- 
ficients, which  are  calculated  using  Equation  28.  The  as- 
sumed forms  of  the  generalized  forces  are  described  in  Sec- 
tion 5.0  on  program  output. 


3.1.3  Kutta  Condition 


The  numerical  treatment  of  the  Kutta  condition  is 
an  essential  part  of  the  calculation  for  lifting  airfoils 
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and  is  worth  discussing  in  some  detail.  Consideration  of 
airfoils  with  camber  or  angle  of  attack  requires  the  cal- 
culation of  the  circulation  (y0)  around  the  airfoil.  This 
is  not  known  apriori  and  must  be  a part  of  the  overall  so- 
lution process.  Application  of  the  Kutta  condition  re- 
quires that  the  circulation  equal  the  jump  in  potential 
across  a "wake"  extending  from  the  airfoil  to  downstream 
infinity.  Incorporation  of  the  wake  in  the  finite  differ- 
ence scheme  involves  differencing  the  4-yy  term  for  points 
near  the  wake  in  such  a manner  as  to  take  into  account  the 
jump  in  <f>.  The  method  for  accomplishing  this,  and  the  it- 
eration method  for  calculating  yo  are  now  discussed. 

The  potentials  just  above  (y  = 0+)  and  below  (y  = 

0 ) the  wake  are  defined  as  cf>+  and  <j>“  respectively,  and 
an  average  potential  at  the  wake  is  defined  in  the  obvi- 
ous manner: 


4>  = *+  (30) 

The  potentials  <f>+  and  4>  satisfy  the  nonlinear  transonic 
potential  equation  and  it  can  easily  be  shown  that  <f>  like- 
wise satisfies  the  same  differential  equation.  Also,  using 
continuity  of  pressure  and  the  potential  equation  it  can  be 
shown  that: 


+ = <p~ 

yy  yy 


(31) 


This  plus  the  fact  that  c|>  = <j>  + a,  where  a is  the  jump 

in  potential  across  the  wake,  provides  a means  for  account- 
ing for  the  Kutta  condition  in  the  finite  difference  ap- 
proximation to  the  potential  equation.  For  grid  points 
on  the  wake  (j  = JW,  i > ITE)  the  average  potential  is 
calculated  and  used  in  the  differencing  of  <j>yy  for  points 
just  above  and  just  below  the  wake.  The  difference  forms 
for  <j>  at  these  points  are  given  by: 


j = JW  + 1 
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yy  i 
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MVi,itr|MViI2j)*i.i+iIM,i  + 2“/  (32) 
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j = JW  (wake) 


♦wli  . = AYVi,j+i  - AY1 

± • J 


j f i + r) 

- AY2j(*i  ' rj  + AYVi,i-i 


(33) 


j = JW  - 1 


yy 


AYl 


1.3 


- V)  - 


(AYl  +AY2  .)<J>  + AY2  <j)  (34) 

J J -WJ  JJ-fJ"-1- 


where  AYl  and  AY2  were  defined  earlier;  and  where 


( x i~  1 ) 

ai  ~ (X(IM)-l)  (Yff“Yte)  + Yte 


(35) 


where  Oi  is  the  jump  in  potential  across  the  wake  at  any  x 
station  and  Yff  and  yte  are  the  circulation  in  the  farfield 
and  at  the  trailing  edge  respectively.  At  any  stage  of  the 
iteration  process,  ai  is  taken  as  the  linearly  interpolated 
value  between  the  current  value  of  yfc  evaluated  as  the 
jump  in  potential  at  the  trailing  edge  (x  = 1) , and  Yff(x  = 

X ( IM) ) . In  the  converged  solution  Yff  = Y.  so  that  S$(x>l)=Y 
is  a constant  along  the  wake  therebyautomaf ically  satisfy- 
ing the  condition  of  pressure  continuity  [<{>]=  0 along  the 
wake . x 


The  differencing  scheme  given  above  is  incorporated 
into  the  process  for  setting  up  the  system  of  algebraic 
equations  for  a column  of  grid  points  with  Xj_  > 1.  The 
scheme  is  symmetric  with  respect  to  points  above  and  below 
the  wake,  and  retains  the  second  order  accuracy  of  the  dif- 
ference approximation  for  fyyy.  <p  on  the  wake  is  calculated 
as  any  other  grid  point  since,  as  noted  earlier,  it  satis- 
fies the  transonic  potential  equation. 


It  remains  only  to  define  a scheme  for  updating  Yff 
based  on  the  calculated  values  of  Y . This  is  done  by 

over-relaxing  on  previous  values  of  y^  . 

te 


ff 


(1) 

Yte 


+ PGFF 


( (NGFF) 
yte 


- Y 


(1) 

te 


(36) 
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where  NGFF  is  the  number  of  grid  sweeps  before  Yff  is  up- 
dated and  PGFF  is  a weighting  factor.  It  is  recommended  that 
NGFF  be  selected  such  that  NGFF  > 1 . 5M  where  M is  the  num- 
ber of  grid  points  between  the  airfoil  and  the  farfield. 

This  guarantees  that  information  from  the  previous  update 
of  the  farfield  has  had  time  to  reach  the  airfoil. 


3.1.4  Farfield  Boundary  Condition 


The  completion  of  the  finite  difference  analog  to 
the  steady  transonic  flow  problem  requires  that  boundary 
conditions  be  fixed  on  the  far  boundary  in  some  manner. 
This  is  accomplished  in  STRANS  in  two  different  ways  de- 
pending upon  the  freestream  Mach  number,  M^.  For  < 1 , 
the  potential  on  all  boundary  points  (I  = 1,  IM  and  J = 1, 
JM)  is  specified  as  a Dirichlet  boundary  condition  and 
for  Mro  > 1 the  small  perturbation  form  of  the  character- 
istic relations  is  used  to  relate  i>  to  <b  . 

ry  Tx 


For  a subsonic  freestream,  analytical  expressions 
are  available  which  are  valid  far  from  the  airfoil.  The 
steady  subsonic  farfield  has  been  studied  most  extensively 
by  Klunker9  and  his  result  is  used  in  STRANS.  The  farfield 
representation  describes  thickness  and  lift  effects  to  lead- 
ing order  and  is  given  by: 
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(37) 


Moo<l,  x2  + Ky2  + °° 


where  t(£)  is  the  airfoil  thickness  distribution.  The 
doublet  strength  due  to  the  nonuniform  flow  is  given  by 
the  integral  of  (4>x)  2 over  the  solution  field.  This  is 
evaluated  in  STRANS  periodically  (every  NGFF  iterations) 
as  the  iterative  solution  proceeds  using  a trapezoidal 
rule  integration  scheme.  Similarly,  in  calculations  for 
airfoils  with  lift,  the  circulation  Yff  is  updated  (also 
every  NGFF  iterations) . Thus  as  the  doublet  strength  and 
airfoil  circulation  are  refined.  Equation  37  is  used  to 
update  <f>ff  on  the  boundaries  of  the  finite  difference  grid. 
Suggestions  for  locating  and  updating  the  farfield  are 
given  in  Section  6.0  below. 
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STRANS  can  also  treat  supersonic  freestream  Mach 
numbers  by  setting  a characteristic  boundary  condition  on 
the  farfield  boundaries.  The  upstream  boundary  is  a com- 
pletely uniform  flow  so  that  (f>°  = 0 is  set  on  i = 1,  j = 1, 
JM.  Because  of  the  backward  differencing  used  for  super- 
sonic grid  points,  no  prescription  is  required  on  the  down- 
stream boundary  as  long  as  all  points  are  supersonic,  as 
they  should  be  in  a converged  solution.  If  during  the  it- 
eration procedure  a grid  point  on  the  boundary  becomes  sub- 
sonic, then  a zero  gradient  (<j>£  = 0)  boundary  condition  is 
used.  On  the  top  and  bottom  of  the  grid,  it  is  assumed 
that  the  disturbance  from  the  body  is  weak  so  that  the 
characteristic  relation  on  the  incoming  characteristic 
applies.  The  first  order  approximation  to  this  relation 
is : 


!Y  (JM) 

(38) 

Y(l) 

This  is  incorporated  in  the  finite  difference  procedure 
in  the  same  manner  as  the  airfoil  boundary  condition. 

That  is,  on  the  top  and  bottom  of  the  grid  (<J>yy)i  j is 
differenced  as  a one-sided  difference  using  Equation  38  in 
the  finite  difference  forms. 


3 . 2 Numerical  Solution  Method  for  First  Order  Unsteady 
Perturbation  Potential 


As  noted  before,  the  boundary  value  problem  for  the 
first  order  unsteady  perturbation  potential  cf> 1 is  of  the 
same  mixed  elliptic/hyperbolic  type  as  the  boundary  value 
problem  for  the  steady  perturbation  potential.  The  govern- 
ing equation  for  <{> 1 is  linear  which  would  suggest  that  many 
of  the  powerful  techniques  of  classical  analysis  could  be 
applied.  The  equation  is  seriously  complicated,  however, 
by  the  fact  that  various  coefficients  are  functions  of  the 
nonlinear  steady  potential,  for  which  the  only  general  solu- 
tion methods  are  numerical.  In  view  of  this  complication 
and  the  success  of  the  mixed  differencing  relaxation  pro- 
cedure for  the  steady  potential  equation  of  similar  type, 
it  was  decided  to  use  this  same  method  in  UTRANS  to  develop 
unsteady  solutions.  The  application  of  the  method  to  the 
unsteady  potential  equation  involves  but  slight  variations 
on  the  procedure  described  in  detail  in  Section  3.1.  These 
variations  are  emphasized  in  the  summary  of  the  solution 
method  given  in  this  section. 
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j = JM 


Figure  5.  Schematic  of  Numerical  Solution  Domain 


Consider  the  boundary  value  problem  for  the  unsteady 
transonic  potential,  with  numerical  solution  domain  indi- 
cated schematically  in  Figure  5.  It  is  recalled  that  the 
governing  equation  is  linear  and,  for  ft  ^ 0,  complex.  Also 
the  equation  is  of  mixed  elliptic/hyperbolic  type  with  the 
local  nature  determined  by  the  steady  perturbation  poten- 
tial <j>°.  Because  of  the  strong  similarity  and  dependence 
of  (J)1  on  <t>°,  a possible  solution  approach  would  be  to  solve 
for  both  simultaneously.  However,  because  of  possible 

different  stability  and  convergence  requirements,  it  was 
decided  to  use  a sequential  approach.  That  is,  since  4> 0 
does  not  depend  on  <f> 1 it  is  solved  independently  by  computer 
program  STRANS,  and  the  resulting  solution  stored  on  mag- 
netic tape.  The  converged  solution  so  obtained,  is  then 
used  in  the  solution  process  for  the  corresponding  <j> 1 . 

This  approach  has  an  addted  benefit  in  that  4>°  need  not  be 
regenerated  for  each  unsteady  boundary  disturbance  of  in- 
terest. 


As  mentioned,  the  numerical  solution  proceeds  in  a 
similar  manner  to  that  for  <f0*  The  local  nature  of  the 
equation  at  each  grid  point  is  determined  by  the  corres- 
ponding value  of  (K-4°)  • • = vi,j  at  the  same  grid  point. 

Then  if  K - <f>£  > 0 (elliptic)  trie  x derivatives  of  <p 1 are 
center  differenced  and  if  K - <f)°  < 0 (hyperbolic)  the 


-20- 


x derivatives  are  backward  differenced.  Near  the  sonic 
line  at  so-called  "parabolic  points,"  (defined  in  Section 
3.1)  the  centered  difference  form  for  Vj_  j and  backward 
difference  form  for  the  x derivatives  of^1  are  used. 

The  finite  difference  approximation  to  the  unsteady  poten- 
tial equation  are  thus  summarized  as  follows: 


. >0,  elliptic 


V.  .(cf)1  )?  . + (4)1  ).  . - [2ift+  ( 0 )E  . ] ( <J> 1 ) ^ . = 0 (39) 
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^ j < 0 , hyperbolic 
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V.  . < 0,  V.  . .>0,  parabolic 
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v.  . (4)1  )H  . + ( 4> 1 ) . . - [2ifM-(<f>°  )H  . ] (<t> 1 ) ? . = 0 (41) 
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where  Vi, j and  Vi,j  and  the  elliptic  (superscript  E)  and 
hyperbolic  (superscript  H)  difference  forms  are  as  defined 
in  Section  3.1  above. 


Using  the  above  forms,  the  finite  difference  equa- 
tions are  set  up  for  each  column  (x  = constant)  in  the 
grid  and  the  resulting  sequence  of  linear  algebraic  equa- 
tions is  solved  by  Gaussian  elimination.  The  equations 
in  this  case  are  linear  so  that  no  iterative  solution  is 
required.  Since  the  ’ s are  in  general  complex  the  above 
process  requires  complex  operations  which  are  straightfor- 
ward on  modern  machines.  After  each  column  is  solved,  it 
is  relaxed  with  a variable  relaxation  factor,  depending  on 
the  local  nature  of  the  solution,  as  before.  This  process 
is  repeated  for  each  column  in  turn,  sweeping  the  grid  from 
left  to  right  until  the  change  in  <J> 1 for  all  grid  points 
during  one  grid  sweep  is  less  than  some  arbitrary  small 
amount  (EPSGRD) . A grid  halving  routine  has  been  imple- 
mented in  UTRANS  in  the  same  manner  as  described  for  STRANS, 
with  similar  improvements  in  the  efficiency  of  the  method. 


-21- 


The  numerical  treatment  of  the  body  boundary  con- 
dition and  Kutta  condition  for  the  unsteady  perturbation 
potential  are  the  same  as  described  earlier  for  <p 0 . The 
only  difference  is  that  all  quantities  such  as  the  cir- 
culation and  jump  in  <p  across  the  wake  are  in  general  com- 
plex. The  iteration  procedure  for  (YjJff  is  the  same  as 
that  given  in  Section  3.1.3  above. 

The  present  version  of  UTRANS  accounts  on  option 
for  two  rigid  body  modes  of  oscillation;  pitch  and  control 
surface  oscillation.  The  specification  of  a non-zero  hinge 
point  indicates  the  desired  mode  is  control  surface  oscil- 
lation. The  respective  airfoil  boundary  conditions  are 
specified  in  subroutine  INITAL  and  are  simply: 


4>y  = -1  on  y = ±0,  xh  < x < 1 (42) 

where  is  the  hinge  point.  Unsteady  components  of  the 
generalized  forces  are  also  calculated  by  a trapezoidal 
rule  integration  of  the  unsteady  pressure  perturbation. 
Complete  definitions  of  the  forces  are  given  in  Section 

5.0. 


The  farfield  boundary  conditions  for  <f> 1 are  defined 
from  an  asymptotic  expression  for  subsonic  freestream  flow 
and  a form  of  the  unsteady  characteristic  relation  for  su- 
personic freestream  flows.  For  Mm  < 1 an  asymptotic  ex- 
pression is  used  to  fix  values  of  °°(J) 1 in  the  Dirichlet  sense 
on  all  four  boundaries  of  the  grid,  thereby  closing  the 
sequence  of  equations  for  <J> 1 at  a column  of  grid  points. 

The  expression  was  derived  in  Reference  3 and  is  given  by: 
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where  I is  given  by: 
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where 


S (R)  = /R2-K£2y2  + R 
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dR 


for  x > 1 


where  k = Q./K , R = £/(x'-£)  2+ky2  and  R1  = k/Tx-1  )'2'+Ky 2 and 
H±  ' is  the  Hankel  function  of  the  second  kind  of  ith  order, 
The  integral  I?  cannot  be  evaluated  in  closed  form  but, 
since  it  is  over  a finite  range,  can  be  integrated  numer- 
ically with  little  difficulty.  It  is  noted  that  the  far- 
field  approximation  depends  on  the  solution  through  the 
unsteady  component  of  the  airfoil  circulation  ( y x ) and  also 
through  the  "airfoil  integral."  The  airfoil  contribution 
to  the  farfield  involves  an  integral  over  the  airfoil 
(°  1 x < 1)  of  the  jump  in  potential  between  the  upper  and 
lower  surface.  Thus  as  the  airfoil  circulation  and  poten- 
tial distribution  is  updated,  the  farfield  is  periodically 
updated  during  the  solution  process,  using  Equation  43. 


Supersonic  freestream  flow  can  also  be  treated  by 
UTRANS  in  much  the  same  manner  as  described  above  for  S TRANS . 
The  appropriate  small  perturbation  approximation  to  the 
characteristic  relation  for  the  unsteady  perturbation  po- 
tential is: 
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This  is  incorporated  into  the  finite  difference  approxi- 
mation to  (j>yy  at  the  top  and  bottom  of  the  grid  using  a 
one  sided  difference  approximation  as  before.  Upstream  and 
downstream  sides  of  the  grid  are  treated  as  described  above. 
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4.0 


PROGRAM  DESCRIPTIONS 


Computer  programs  STRANS  and  UTRANS,  used  in  con- 
junction, implement  the  theory  and  numerical  solution  pro- 
cedure for  unsteady  transonic  flow  described  in  the  previ- 
ous two  sections.  As  described  above,  the  boundary  value 
problems  for  the  steady  perturbation  potential  (Equation  10) 
and  the  unsteady  perturbation  potential  (Equation  11)  are 
solved  in  STRANS  and  UTRANS  respectively  using  a finite 
difference  relaxation  procedure.  The  use  and  manipulation 
of  magnetic  tapes  forms  an  integral  part  of  the  operation 
of  each  program  as  well  as  serving  as  the  necessary  "data 
link"  between  the  two  programs.  As  a result  the  user  is 
assumed  to  have  some  familiarity  with  the  use  of  tapes  and 
their  manipulation  with  control  cards.  The  reading  and 
writing  of  data  files  on  magnetic  tape  is  described  in  the 
next  section  and  motivated  in  Section  6.0.  In  this  section, 
the  logical  flow  of  the  STRANS  and  UTRANS  programs  is  des- 
cribed and  a brief  summary  of  each  subroutine  is  presented. 
Both  programs  are  quite  similar  in  logical  approach  and 
operation,  so  that  they  are  described  together.  Differ- 
ences between  the  programs  are  highlighted  with  appropri- 
ate comments  as  needed. 

The  logical  flow  of  the  STRANS  program  is  shown 
schematically  in  Figure  6 below.  The  general  logic  of  the 
UTRANS  program  is  almost  identical  with  minor  exceptions 
noted  in  the  description  below.  The  calculation  is  begun 
by  reading  card  input  and,  if  a restart  is  being  performed, 
a tape  dump.  In  UTRANS  the  tape  dump  of  the  steady  solu- 
tion being  perturbed  is  also  read.  All  finite  difference 
coefficients  and  airfoil  boundary  conditions  are  initial- 
ized in  a call  to  INITAL  and  subsonic  farfield  quantities 
are  initialized  in  a call  to  FARFLD.  If  a restart  is  not 
being  performed,  initial  values  for  <j>  at  all  grid  points 
are  determined  by  the  linearized  subsonic  or  supersonic 
solution.  The  computational  cycle  is  executed  by  setting 
up  the  tridiagonal  equations  for  a column  of  grid  points 
using  the  mixed  differencing  finite  difference  equations. 

The  equations  are  solved  iteratively  in  STRANS  and  in  one 
pass  in  UTRANS,  by  Gaussian  ellimination  in  a call  to  TRI . 
Each  column  is  solved  and  relaxed  in  turn  proceding  through 
the  grid  from  left  to  right.  The  grid  is  swept  iteratively 
in  this  manner  until  the  change  in  <J>  for  all  grid  points  is 
less  than  EPSGRD(l).  A call  to  PRINT  prints  out  the  air- 
foil pressure  coefficients  every  NPRINT  iterations,  and 
the  farfield  is  updated  every  NGFF  iterations,  in  FARFLD. 
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STRANS 


• Driver  Routine 

• Continuous  Commentary 

• Calls  to  Other  Routines 
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Figure  6, 


Logical  Flow  of  STRANS  Program 
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Figure  6.  Logical  Flow  of  STRAUS  Program  (Continued) 
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Figure  6 


. Logical  Flow  of  STRANS  Program  (Continued) 
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When  the  converged  solution  is  obtained,  a tape  dump 
of  all  relevant  input  and  calculated  quantities  is  performed 
and  a call  to  FPRINT  calculates  and  prints  out  the  airfoil 
pressure  and  force  coefficients.  Various  diagnostic  prints 
are  also  performed  in  STRANS  and  UTRANS  after  every  grid 
iteration,  when  the  farfield  is  updated,  when  the  grid  is 
refined  and  when  a tape  dump  is  performed.  Finally,  if  de- 
sired, the  grid  is  refined  by  doubling  all  grid  lines  in  a 
call  to  subroutine  D0UBLE , and  the  computational  cycle  is 
begun  anew  until  the  new  solution  has  again  converged  or 
until  the  maximum  number  of  iterations  (NGRID)  has  been  ex- 
ceeded. In  either  case,  a final  tape  dump  and  final  print 
are  executed. 

A summary  of  each  subroutine  is  now  presented. 


STRANS/UTRANS 


These  are  the  driver  routines  for  the  respective  pro- 
grams. The  logical  flow  of  the  mixed  differencing  relaxa- 
tion procedure  as  just  described  is  controlled  by  these 
routines  and  all  operations  including  input,  initialization, 
finite  difference  solution  and  output  are  performed  either 
internally  or  by  calls  to  the  various  subroutines  described 
below. 


D0UBLE 


This  routine  refines  the  grid  spacing  by  doubling 
the  number  of  grid  lines.  New  grid  points  are  added  in 
both  the  X and  Y direction  by  placing  grid  lines  midway  be- 
tween the  existing  ones.  Initial  values  for  <f>  at  the  new 
grid  points  are  determined  by  a linear  interpolation.  Calls 
are  made  to  INITAL  and  FARFLD,  to  reevaluate  finite  dif- 
ference coefficients  and  regenerate  farfield  values  for  4> 
respectively . 


FARFLD 


The  subsonic  farfield  is  calculated  and  updated  in 
this  routine  using  the  asymptotic  solutions  for  the  steady 
(Equation  37)  or  unsteady  (Equation  43)  perturbation  poten- 
tials. In  the  STRANS  version,  the  doublet  strength  (D0UBLT ) 
is  also  updated  here  by  performing  a trapezoidal  rule  inte- 
gration over  the  finite  difference  grid  of  the  expression: 
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D0UBLT  = D0UB  + 


(<J>£)  2dxdy 


FL  (in  STRANS  only) 

This  is  a function  statement  which  contains  the  air- 
foil lower  surface  shape  function.  This  function  appears 
in  the  supersonic  small  disturbance  solution  which  is  used 
in  STRANS  as  an  initial  guess  for  <j>°  below  the  airfoil  for 
supersonic  freestream  flows. 


FLP  (in  STRANS  only) 

This  is  a function  statement  which  contains  the  air- 
foil lower  surface  slope  distribution  used  in  the  linearized 
tangency  boundary  condition.  This  function  is  called  from 
subroutine  INITAL  and  its  value  at  each  grid  point  on  the 
lower  surface  of  the  airfoil  is  stored  in  the  FPL  array. 


FPRINT 


This  routine  produces  the  final  print  and  is  called 
when  the  solution  has  converged  to  the  desired  accuracy  or 
when  the  problem  is  terminated  for  reaching  the  maximum  num- 
ber of  grid  iterations  allowed  (NGRID) . The  unsealed  pres- 
sure coefficients  above  and  below  the  airfoil  and  the  air- 
foil force  coefficients  are  also  calculated  and  printed  out 
in  this  routine. 


FU  (in  STRANS  only) 

This  is  a function  statement  which  contains  the  air- 
foil upper  surface  shape  function.  This  function  appears 
in  the  supersonic  small  disturbance  solution  which  is  used 
in  STRANS  as  an  initial  guess  for  cp°  above  the  airfoil  for 
supersonic  freestream  flows. 


FUP  (in  STRANS  only) 

This  is  a function  statement  which  contains  the  air- 
foil upper  surface  slope  distribution  used  in  the  linearized 
tangency  boundary  condition.  This  function  is  called  from 
subroutine  INITAL  and  its  value  at  each  grid  point  on  the 
upper  surface  of  the  airfoil  is  stored  in  the  FPU  array. 
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The  doublet  strength  due  to  airfoil  thickness  ( D0UB ) must 
also  be  given  in  this  subroutine.  This  quantity  is  de- 
fined by  an  integral  of  the  airfoil  thickness  distribution 
function  (normalized  to  airfoil  thickness) : 


This  routine  performs  the  extrapolation  to  update 
farfield  circulation  (GAMFF ) . 


HANKEL  (in  UTRANS  only) 

This  subroutine  calculates  the  first  and  second 
order  Hankel  functions  of  the  second  kind  for  real  argu- 
ments using  curvefits  developed  in  Reference  10.  These 
functions  are  used  in  the  calculation  of  the  farfield 
for  the  unsteady  perturbation  potential. 


The  finite  difference  coefficients  AX1 , AX2 , BX1, 
BX2 , CX,  AYl  and  AY2  and  AX(DX)  and  AY(DY)  are  computed 
in  this  subroutine.  The  airfoil  boundary  conditions  FPU 
and  FPL  are  also  set  here,  using  functions  FUP  and  FLP 
respectively. 


This  routine  computes  and  prints  the  scaled  pres- 


sure coefficients  above  and  below  the  airfoil  every  NPRINT 
grid  iterations. 


This  routine  solves  a system  of  tridiagonal  equa- 
tions using  Gaussian  elimination. 


o 


GAMFUN 


INITAL 


PRINT 


TRI 
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5.0 


INPUT  AND  OUTPUT 


A description  of  the  input  required  to  run  STRANS 
and  UTRANS,  and  the  output  associated  with  each  program 
are  presented  in  this  section.  All  card  input  is  entered 
using  the  standard  CDC  NAMELIST  package  with  the  exception 
of  a title  card. 

Figure  7 below  is  presented  to  summarize  the  defi- 
nitions of  various  important  parameters  used  in  STRANS  and 
UTRANS.  It  is  reiterated  that  STRANS  calculates  the  steady 
transonic  flow  about  a general  thin  airfoil  with  mean  angle 
of  attack  a0  and  with  a control  surface  with  hinge  point 
and  mean  flap  angle  a|.  UTRANS  calculates  the  unsteady 
perturbation  to  this  steady  flow  consisting  of  an  harmonic 
oscillation  in  angle  of  attack  about  the  nose  or  of  a con- 
trol surface  oscillation.  The  output  parameters  of  great- 
est interest  are  the  airfoil  pressure  and  force  coefficients 
which  are  printed  out  once  the  final  converged  solution  is 
obtained.  The  definitions  of  the  scaled  and  unsealed  forms 
of  the  pressure  coefficients  are  shown  in  Figure  7.  Also 
defined  in  the  figure  are  the  relations  used  in  the  programs 
for  the  steady  and  unsteady  components  of  the  generalized 
lift,  moment  about  x = 0 and  hinge  moment  coefficients.  It 
is  noted  that  counterclockwise  moments  are  defined  positive 
and  that  the  unsteady  perturbations  to  all  coefficients  are 
given  per  unit  angle  of  oscillation.  The  unsteady  pertur- 
bations to  pressure  and  force  coefficients  are  in  general 
imaginary  so  that  the  real  and  imaginary  components  are 
printed  out  in  the  final  print  of  UTRANS. 
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Figure  7,  Definition  of  Pressure  and  Force  Coefficients 


5.1 


STRANS  Input 


The  input  for  STRANS  is  now  considered  in  three  sets. 
Recommended  and/or  typical  values  for  some  of  the  input  var- 
iables which  control  the  numerical  scheme,  appear  in  paren- 
theses. Also  presented  at  the  end  of  this  section  is  a 
description  of  the  restart  capability  which  requires  input 
from  a magnetic  tape  dump  of  a previous  calculation. 

First  Set 


BCD  title  card  containing  any  information  in  columns 
1 through  80  (Format  8A10) . This  can  be  used  to  define  the 
case  being  run  and  is  printed  out  on  the  last  page  of  output 
which  presents  the  final  converged  results. 

Second  Set 


The  second  set  of  data  is  read  in  under  NAMELIST 
name  $C0NTRL.  The  single  variable  read  defines  the  use 
of  the  restart  option.  Some  comments  concerning  the  me- 
chanics of  the  use  of  this  option  are  given  at  the  end  of 
the  section. 


NAME 


DESCRIPTION 


ITAPE  This  is  a flag  for  using  a restart  tape. 

ITAPE  = 0 means  the  problem  is  being  started 
from  scratch  (iteration  0)  using  an  initial 
guess  defined  in  STRANS.  ITAPE  = 1 means 
the  problem  is  being  restarted  from  a pre- 
vious run  which  is  to  be  read  from  a dump 
tape. 

Third  Set 


The  third  set  of  data  is  read  in  under  NAMELIST 
name  $IN,  and  includes  all  of  the  variables  required  to 
define  a problem,  and  control  the  numerical  iteration 
procedure . 
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name  description 


X 

An  array  containing  the  streamwise  grid 
coordinates;  IM  of  them 

Y 

An  array  containing  the  normal  grid  coord- 
inates; JM  of  them 

IM 

Number  of  grid  points  in  the  streamwise 
direction  (maximum  of  100) 

JM 

Number  of  grid  points  in  the  normal  direc- 
tion (maximum  of  100) 

ILE 

I location  of  airfoil  leading  edge  (X(ILE)) 

ITE 

I location  of  airfoil  trailing  edge  (X(ITE)) 

JW 

J location  of  airfoil  (Y(JW)) 

M8 

Freestream  Mach  number 

GAM 

y,  ratio  of  specific  heats 

DEL 

Airfoil  thickness  ratio  in  percent 

ALPHA 

Airfoil  angle  of  attack  in  radians 

ALPHAF 

Flap  angle  in  radians 

XH 

X location  of  the  hinge  point  on  the  airfoil 
(0  £ XH  < 1) 

GAMFF 

Initial  guess  for  the  airfoil  circulation; 
to  be  used  in  the  initialization  of  the 
farf ield 

NGFF 

Every  NGFF  grid  iterations  the  farf ield  cir- 
culation and  doublet  strength  are  updated. 
This  also  causes  the  farf ield  to  be  updated. 
(^10) 

PGFF 

Relaxation  parameter  used  in  the  iteration 
for  the  airfoil  circulation  in  the  farfield 
(^1.5) 

0MEGAH 

Relaxation  parameter  for  hyperbolic  grid 
points  (^.75) 
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NAME  DESCRIPTION 


0MEGAE 

Relaxation  parameter  for  elliptic  grid  points 
(VL.7) 

0MEGAP 

Relaxation  parameter  for  parabolic  grid 
points  (^.75) 

EPSC0L 

Convergence  criteria  for  column  solution. 

The  change  in  |<J>°|  during  a column  iteration 
at  every  point  in  the  column  must  be  less 
than  EPSC0L  for  convergence  to  occur 
(^5  x 10-5) 

NC0L 

Maximum  number  of  column  iterations  allowed. 
Note  that  if  NC0L  iterations  is  reached 
without  convergence,  a printout  of  the  degree 
of  convergence  is  given  and  the  calculations 
proceed  as  if  convergence  had  occurred  (^10) 

EPSGRD 

An  array  containing  criteria  to  control 
grid  convergence.  The  change  in  |<j>°|  at 
every  grid  point  during  one  grid  sweep 
must  be  less  than  EPSGRD (1)  for  conver- 
gence to  occur.  The  number  of  grid  lines 
is  then  doubled  if  KEPS  > 1 and  EPSGRD (2) 
is  used  to  determine  grid  convergence,  and 
so  on.  KEPS  values  of  EPSGRD  are  input. 

KEPS 

The  number  of  EPSGRD ' s input  (maximum  of  3). 
The  grid  is  refined  KEPS-1  times. 

NGRID 

Maximum  number  of  grid  iterations  allowed. 
When  the  number  of  grid  iterations  equals 
NGRID  the  calculation  is  terminated. 

NDUMP 

Binary  tape  dump  frequency.  Every  NDUMP 
grid  iterations  current  values  of  all  vari- 
ables will  be  dumped  on  tape.  Note  that  a 
tape  dump  occurs  automatically  whenever  the 
grid  converges  or  the  number  of  grid  itera- 
tions equals  NGRID  (set  equal  to  large  num- 
ber if  a dump  of  only  the  final  iteration  is 
desired) . 

NPRINT 

Every  NPRINT  grid  iterations  the  scaled 
pressure  coefficient  above  and  below  the 
airfoil  is  printed. 
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NAME 


DESCRIPTION 


IK  Setting  IK  = 1,  flags  the  use  of  a previous 

solution  as  an  initial  guess  for  the  current 
problem  where  the  Mach  number,  airfoil  thick- 
ness or  shape  and/or  angle  of  attack  may  be 
different . 


The  input  data  listed  above  are  necessary  to  initi- 
ate a calculation  for  which  no  initial  guess  is  available. 
Most  calculations,  however,  are  performed  as  restarts  using 
data  which  has  been  stored  as  binary  files  on  the  restart 
tape  according  to  the  format  described  in  Section  5.2. 

This  use  of  the  restart  capability  is  an  inherent  aspect 
of  the  recommended  computational  procedure  as  is  discussed 
in  some  detail  in  Section  6.0.  Some  brief  comments  des- 
cribing the  initiation  of  a calculation  using  the  restart 
capability  are  pertinent  at  this  juncture. 

It  is  noted  that  the  restart  or  dump  tape  ( TAPE7 ) 
may  be  manipulated  in  any  way  desired  using  the  appropri- 
ate control  cards.  In  general  the  tape  will  contain  data 
from  many  runs,  stored  as  individual  binary  files.  For 
restarting  the  STRANS  program,  the  desired  file  from  the 
restart  tape  (TAPE7 ) is  copied  to  a disc  file  (TAPE8 ) . 

The  user  is  reminded  to  rewind  TAPE 8 . TAPE 7 is  then  po- 
sitioned at  the  end  of  the  last  file  on  the  tape  so  that 
new  dumps  can  be  written  by  the  program  without  losing  any 
of  the  old  data.  The  first  two  sets  of  data  are  then  in- 
put with  ITAPE=1.  In  the  third  set  of  data  the  following 
control  variables  are  needed  as  input: 

0MEGAH , 0MEGAE , 0MEGAP , EPSC0L , EPSGRD , NDUMP , 

NC0L , NGRID , NGFF , PGFF , KEPS  and  NPRINT . 

The  remaining  input  variables  are  stored  on  the  restart 
tape  and  need  not  be  input  unless  the  restart  option  is 
being  used  to  run  a new  case.  If  a new  case  is  being  run, 

IK  must  be  set  to  1 which  allows  the  Mach  number,  airfoil 
thickness,  airfoil  angle  of  attack,  flap  angle  and/or  the 
hinge  position  (M8,  DEL,  ALPHA,  ALPHAF  and  XH)  to  be  changed. 
If  the  airfoil  angle  of  attack  and/or  the  flap  angle  are 
changed  a new  guess  for  the  farfield  circulation  (GAMFF ) 
can  and  should  be  made. 
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5.2  STRANS  Output 


The  output  from  STRANS  consists  of  three  parts: 

(i)  a continuous  commentary  which  describes  the  progress 
of  the  iterative  solution  procedure,  (ii)  a final  print 
summarizing  results  of  interest  from  the  final  converged 
solution,  and  (iii)  a binary  tape  dump  of  all  pertinent 
input  and  calculated  parameters. 

The  continuous  commentary  consists  of  various  print 
statements  executed  in  the  main  program  STRANS  or  the  sub- 
routine PRINT  which  describe  the  current  state  of  the  solu- 
tion as  well  as  the  occurrence  of  various  "milestones"  in 
the  iteration  process.  The  only  print  that  occurs  every 
iteration  is  the  value  of  the  maximum  change  in  <J>°  through- 
out the  grid  during  one  grid  iteration.  When  a column  it- 
eration fails  to  converge,  a print  occurs  which  defines  the 
degree  of  column  convergence  and  the  j location  of  the  most 
poorly  converged  point.  Every  NGFF  iterations,  the  subsonic 
farfield  is  updated  and  the  new  values  of  farfield  circulation 
(GAMFF ) , airfoil  circulation  (GAMTE)  and  nonlinear  doublet 
strength  (D0UBLT)  are  printed.  The  user  can  examine  the  ef- 
fect of  degree  of  convergence  on  the  solution  by  specifying 
a print  of  the  scaled  pressure  coefficients  on  the  upper 
and  lower  airfoil  surfaces  every  NPRINT  iterations.  Finally, 
a descriptive  print  occurs  at  certain  milestone  points  such 
as  the  doubling  of  grid  lines,  the  occurrence  of  a binary 
tape  dump  and  solution  convergence. 

The  final  print  is  executed  in  subroutine  FPRINT 
when  the  solution  has  converged  to  the  desired  accuracy 
or  when  the  number  of  grid  iterations  equals  NGRID.  The 
print  is  self-explanatory  and  includes  the  input  parameters 
which  define  the  problem  and  various  calculated  quantities 
of  interest.  The  calculated  quantities  are  of  course  based 
on  the  final  converged  solution.  The  steady  force  coeffi- 
cients, as  defined  in  Figure  7,  are  printed  out  as  well  as 
the  upper  and  lower  surface  pressure  coefficients  with  cor- 
responding airfoil  coordinate. 

The  most  important  form  of  STRANS  output  is  the  bi- 
nary tape  dump  of  all  input  parameters  defining  the  prob- 
lem and  of  the  most  recent  values  of  <J> 0 at  all  grid  points. 

A tape  dump  occurs  automatically  if  the  solution  has  con- 
verged to  the  desired  accuracy  or  if  the  number  of  grid 
iterations  equals  NGRID.  The  user  may  also  specify  that 
such  a dump  occur  every  NDUMP  grid  iterations.  The  tape 
so  generated,  not  only  forms  a permanent  record  of  the 
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results  of  a calculation  for  possible  future  editing  and 
examination  but  also  forms  a necessary  part  of  the  compu- 
tational procedure.  Most  important  is  its  use  as  required 
input  for  a UTRANS  calculation.  However,  it  may  also  be 
used  to  restart  the  calculation  to  refine  accuracy  or  con- 
vergence or  be  used  as  the  initial  guess  for  <j)°  throughout 
the  grid  for  a similar  calculation,  as  described  in  Sec- 
tion 5.1. 

The  format  used  for  writing  and  reading  the  binary 
tape  is  given  in  the  following  FORTRAN  statements: 


WRITE 

(7) 

NITERG,  IM,  IMl,  JM,  JMl , JW,  JWPl, 
JWM1,  ITE,  ILE , K,  DEL,  ALPHA, 
GAMTEl , GAMFF,  NDB , XH,  M8 , GAM, 
DYBUl , DYBU2 , DYBL1 , DYBL2 , D0UB , 
D0UBLT , ALPHAF 

WRITE 

(7) 

(X(I) ,1=1,  IM) 

WRITE 

(7) 

( Y ( I ) ,1=1,  JM) 

WRITE 

(7) 

(FPU (I),  FPL (I),  PHIUB(I),  I = ILE, 
ITE) 

WRITE 

(7) 

(AX1(I),  AX2  (I ) , BX1(I),  BX2  ( I ) , 
CX(I) , DX ( I ) , I = 1,  IMl) 

WRITE 

(7) 

(AY1(I),  AY2 (I ) , DY  ( I ) , I = l,JMl) 

WRITE 

(7) 

(PHI (I) , I = 1 , L) 

END  FILE  7 

where  IMl  = IM-1,  JM1  = JM-1,  L = IM*JM 


Any  information  may  be  retrieved  from  the  tape  by  using 
the  appropriate  READ  statements  as  is  done  in  the  restart 
option  described  above. 


5 . 3 UTRANS  Input 


The  input  for  UTRANS  consists  of  normal  card  input 
plus  input  from  a binary  file  which  contains  the  steady 
solution  generated  by  an  STRANS  run.  UTRANS  also  has  a 
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restart  capability  which  is  implemented  in  exactly  the  same 
manner  as  previously  described  in  Section  5.1  for  STRANS 
and  elaborated  upon  at  the  end  of  this  section.  The  re- 
quired input  is  now  described  and  some  comments  are  pre- 
sented at  the  end  of  this  section  pertaining  to  the  tape 
read  of  the  steady  solution.  As  before,  the  card  input 
is  described  in  three  sets. 

First  Set 


BCD  title  card  containing  any  information  in  columns 
1 through  80  (Format  8A10) . 

Second  Set 

The  second  set  of  data  is  read  in  under  NAMELIST 
name  $C0NTRL . 


NAME  DESCRIPTION 


ITAPE  This  is  a flag  for  using  a restart  tape. 

ITAPE=0  means  the  problem  is  being  started 
from  scratch  (iteration  0) , ITAPE=1  means 
the  problem  is  being  restarted  from  a pre- 
vious run  using  the  restart  tape.  Note  that 
a tape  is  also  used  for  the  input  of  steady 
results  independent  of  the  value  of  ITAPE. 

Third  Set 


The  third  set  of  data  is  read  in  under  NAMELIST 
name  $IN. 


NAME 


DESCRIPTION 


X 


An  array  containing  the  streamwise  grid 
coordinates;  IM  of  them. 


Y An  array  containing  the  normal  grid  coordi- 

nates; JM  of  them. 

IM  Number  of  grid  points  in  the  streamwise  di- 

rection (maximum  of  100). 
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NAME  DESCRIPTION 


JM 

Number  of  grid  points  in  the  normal  direc- 
tion (maximum  of  100) . 

ILE 

I location  of  airfoil  leading  edge  (X(ILE)). 

ITE 

I location  of  airfoil  trailing  edge  (X(ITE)). 

JW 

J location  of  airfoil  (Y (JW) ) . 

SMALLK 

Reduced  frequency  based  on  chord  = ooc/U. 

XH 

X location  of  the  hinge  point  on  the  air- 
foil. Set  XH=0  for  oscillation  in  pitch 
(0  £ XH  < 1) . 

GAMFF 

Initial  guess  for  the  airfoil  circulation 
used  in  the  initialization  of  the  farfield. 
Note  that  GAMFF  is  a complex  number. 

NGFF 

Every  NGFF  grid  iterations  the  airfoil  cir- 
culation in  the  farfield  is  updated.  This 
also  causes  the  farfield  to  be  updated  (^10) . 

PGFF 

Relaxation  parameter  used  in  the  iteration 
for  the  airfoil  circulation  in  the  farfield 
(VL.5)  . 

0MEGAH 

Relaxation  parameter  for  hyperbolic  grid 
points  (^ . 75) . 

0MEGAE 

Relaxation  parameter  for  elliptic  grid 
points  (^1 . 7 ) . 

0MEGAP 

Relaxation  parameter  for  parabolic  grid 
points  ('v.  75)  . 

EPSGRD 

An  array  containing  criteria  to  control  grid 
convergence.  The  change  in  | <f> 1 | at  every  grid 
point  during  one  grid  sweep  must  be  less  than 
EPSGRD (1)  for  convergence  to  occur.  The  num- 
ber of  grid  lines  is  then  doubled  if  KEPS  > 1 
and  EPSGRD  (2)  is  used  to  determine  grid  conver- 
gence and  so  on.  KEPS  values  of  EPSGRD  are 
input. 

KEPS 

The  number  of  EPSGRD ' s input  (maximum  of  3). 
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NAME 


DESCRIPTION 


NGRID 


NDUMP 


NPRINT 


IK 


Maximum  number  of  grid  iterations  allowed. 
When  the  number  of  grid  iterations  equals 
NGRID  the  calculation  is  terminated. 

Binary  tape  dump  frequency.  Every  NDUMP 
grid  iterations  current  values  of  all  vari- 
ables will  be  dumped  on  tape.  Note  that  a 
tape  dump  occurs  automatically  whenever 
the  grid  converges  or  the  number  of  grid 
iterations  equals  NGRID.  (Set  equal  to 
large  number  if  a dump  of  only  the  final 
iteration  is  desired.) 

Every  NPRINT  grid  iterations  the  scaled 
upper  and  lower  surface  pressure  coefficient 
per  unit  angle  of  oscillation  is  printed. 

Setting  IK=1  allows  the  user  to  use  a previ- 
ous solution  as  an  initial  guess  for  the  cur- 
rent problem  where  the  reduced  frequency 
and/or  mode  of  oscillation  is  different. 


It  is  recalled,  that  the  solution  for  the  unsteady 
perturbation  potential,  implemented  in  UTRANS,  requires 
the  solution  of  the  steady  potential,  generated  by  STRANS. 
This  is  accomplished  by  reading  the  appropriate  file  on  a 
dump  tape  generated  by  STRANS,  in  much  the  same  way  as  is 
done  in  the  restart  option.  It  is  instructive  to  briefly 
describe  the  UTRANS  restart  including  the  tape  read  of  the 
steady  solution. 

Restarting  the  UTRANS  program  is  only  slightly  more 
complicated  than  STRANS.  In  this  case,  two  tape  dumps  or 
files  are  required.  First  the  file  containing  the  desired 
steady  tape  dump  is  copied  from  TAPE7  to  a disc  file,  TAPE8. 
Next  the  file  containing  the  desired  unsteady  tape  dump  is 
copied  from  TAPE 7 to  a disc  file,  TAPE9 , and  TAPE 8 and  TAPE 9 
are  rewound  for  reading  by  the  program.  TAPE7  is  then  po- 
sitioned at  the  end  of  the  last  file  on  the  tape  in  prepara- 
tion for  accepting  a new  tape  dump.  The  first  two  sets  of 
data  are  input  as  before  (be  sure  to  set  ITAPE=1) . In  the 
third  set  of  data  the  following  variables  are  necessary: 

0MEGAH , 0MEGAE , 0MEGAP , EPSGRD , NDUMP,  NGRID, 

NGFF , PGFF , KEPS  and  NPRINT. 
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Again  there  is  an  option  (IK=1)  which  allows  the  user  to 
change  the  reduced  frequency  and/or  the  mode  of  oscilla- 
tion (SMALLK  and  XH) . In  either  case  a new  guess  for  the 
farfield  circulation  (GAMFF ) should  be  made. 


5.4  UTRANS  Output 


The  output  from  UTRANS  is  very  similar  to  that  of 
STRANS  and  includes  a continuous  commentary,  final  print 
and  binary  tape  dump. 

The  printed  output  is  of  the  form  described  above 
for  STRANS . The  only  difference  is  that  the  field  varia- 
bles in  UTRANS  are  complex  so  that  the  real  and  imaginary 
parts  are  printed  out  in  that  order.  The  descriptive  prints 
are  all  the  same  with  the  deletion  of  the  unneeded  comment 
on  column  convergence.  The  final  print  is  executed  in  sub- 
routine FPRINT  when  the  solution  has  converged  or  has 
reached  the  maximum  number  of  iterations  desired  by  the 
user  (NGRID) . The  print  includes  all  important  input  var- 
iables which  define  both  the  steady  solution  being  perturbed 
and  the  unsteady  solution  being  generated.  Also,  various 
calculated  quantities,  based  on  the  final  converged  solu- 
tion, are  printed.  These  include  the  real  and  imaginary 
parts  of  the  unsteady  contribution  (per  unit  angle  of  os- 
cillation) to  the  aerodynamic  force  coefficients  as  defined 
in  Figure  7 above.  Also  unsteady  contributions  to  the  up- 
per and  lower  surface  pressure  coefficients  (per  unit  angle 
of  oscillation)  are  printed  for  every  computational  point 
on  the  airfoil.  It  is  again  noted  that  these  are  complex 
so  that  the  real  and  imaginary  parts  are  printed  out  in 
order . 


The  other  form  of  UTRANS  output  is  the  binary  tape 
dump  of  all  input  parameters  and  the  most  recent  values  of 
(Reef)1,  Im<f> 1 ) at  all  grid  points.  As  before  the  tape  dump 
occurs  automatically  at  normal  program  termination  or  at 
the  users  discretion  every  NDUMP  iterations.  The  format 
used  for  writing  and  reading  the  binary  tape  is  given  in 
the  following  FORTRAN  statements: 
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WRITE 

(7) 

NITERG,  IM,  IMl,  JM,  JMl,  JW, 

JWP1,  JWMl , ILE,  ITE,  GAMTE1, 

GAMFF , 0MEG , SMALLK,  DYBUl , DYBU2 , 
DYBLl , DYBL2,  ND0UB,  XH 

WRITE 

(7) 

(X  (I)  ,1  = 1,  IM) 

WRITE 

(7) 

( Y ( I ) ,1  = 1,  JM) 

WRITE 

(7) 

(FPU (I ) , FPL ( I ) , PHIUB(I),  I = 
ILE,  ITE) 

WRITE 

(7) 

(AX1(I),  AX 2(1),  BX1(I),  BX2 ( I ) , 
CX(I) , DX ( I ) , I = 1,  IMl) 

WRITE 

(7) 

( AY 1(1),  AY 2 ( I ) , DY ( I ) , I = 1,  JMl) 

WRITE 

(7) 

(PHI  (I)  , I = 1,  L) 

WRITE 

(7) 

(C0NAFF ( I ) , PHIAFF ( I ) , 1=1,  NPT) 

END  FILE  7 

where  IMl  = IM-1,  JMl  = JM-1,  L = NPT  = 

2 * ( IM-3 ) +3*JM 


It  is  again  noted  that  the  above  formats  can  be  used  to 
retrieve  data  from  the  dump  tape  using  the  appropriate 
READ  statements. 
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6.0 


PROGRAM  USAGE 


A summary  of  the  basic  information  needed  for  using 
the  STRANS  and  UTRANS  programs  and  some  hints  on  setting 
up  and  executing  a calculation  are  presented  in  this  sec- 
tion. The  comments  given  here  are  incorporated  in  the  sam- 
ple cases,  presented  in  the  next  section. 

By  way  of  introduction,  it  should  be  noted  that  the 
programs  are  not  now  and  may  never  be  of  a degree  of  so- 
phistication for  which  they  may  be  treated  as  "black  boxes" 
for  the  calculation  of  transonic  aerodynamic  coefficients. 
Finite  difference  numerical  schemes  for  nonlinear  problems 
such  as  transonic  flow  are  not  yet  sufficiently  advanced 
to  permit  the  calculation  of  general  cases  without  the 
care,  attention  and  often  agonizing  of  the  user.  Consid- 
erable judgment  is  required  of  the  user  in  exercising  the 
programs,  and  this  can  only  come  through  experience  with 
similar  calculations.  Section  6.1  is  included  with  the 
hope  of  providing  some  guidance  concerning  various  impor- 
tant computational  details.  In  addition.  Section  6.2  pro- 
vides a summary  of  recommended  procedures  for  exercising 
the  programs  to  provide  the  first  time  user  with  some 
baseline  experience  level  based  on  calculations  performed 
to  date  with  the  programs. 


6 . 1 Computational  Details 


In  this  section  the  general  structure  of  the  com- 
puter programs  is  first  discussed  followed  by  comments  on 
a variety  of  details  for  setting  up  and  running  a calcu- 
lation . 


In  their  present  configuration,  both  STRANS  and  UTRANS 
allow  a maximum  of  10000  computational  grid  points  and  the 
number  of  grid  lines  in  the  streamwise  and  normal  directions 
must  each  be  less  than  or  equal  to  100.  In  this  configura- 
tion, STRANS  requires  638K  words  to  load  and  51eK  words  to 
execute  while  UTRANS  requires  144sK  words  to  load  and  1328K 
words  to  execute.  Greater  storage  is  required  by  UTRANS 
because  the  unsteady  potential  is  complex  so  that  Re^1  and 
Im<J) 1 must  be  stored  in  addition  to  storage  required  for  the 
<J)°  solution  being  perturbed.  The  programs  have  been  exer- 
cised to  date  on  CDC  6600  and  7600  computers.  Based  on 
calculations  performed  to  date  on  the  CDC  7600,  using  the 
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extended  FORTRAN  compiler  (FTN)  with  full  optimization 
(0PT=2 ) , typical  execution  times  are  4.8  x 10“ 5 CPU  sec/ 
grid  point/iteration  for  STRANS  and  6.6  x 10“ 5 CPU  sec/ 
grid  point/iteration  for  UTRANS.  Execution  times  on  the 
CDC  6600  using  the  same  compiler  are  generally  a factor 
of  4 greater.  The  execution  times  given  here  are  useful 
for  estimating  computer  time  requirements  for  a calculation. 

Specification  of  Airfoil  Shape:  All  quantities 

required  to  perform  an  STRANS  calculation  are  input  to 
the  program  via  cards  or  magnetic  tape  (restart)  with  the 
exception  of  the  airfoil  shape.  As  indicated  in  Section 
4.0,  the  airfoil  upper  and  lower  surface  shape  functions 
(normalized  to  the  airfoil  thickness)  are  given  in  function 
subprograms  FU  and(FL  respectively  and  the  upper  and  lower 
surface  slopes  (f “ ' (x)  and  f^'U))  are  given  in  FUP  and  FLP 
respectively.  FU  and  FL  are  needed  only  for  initial  calcula- 
tions with  a supersonic  freestream  and  are  used  in  the 
initial  guess  for  values  of  <j> 0 in  the  grid.  f ° ' and  f | ' and 
the  doublet  strength  due  to  thickness,  D0UB,  need  to  be 
specified  for  both  subsonic  and  supersonic  freestream  cases 
as  they  are  used  in  the  airfoil  boundary  condition.  The 
functions  presently  programmed  in  STRANS  are  for  a symmet- 
ric circular  arc  airfoil.  Thus  to  run  calculations  for 
another  airfoil,  the  respective  statements  need  to  be 
changed  in  the  function  subprograms.  The  specification  of 
airfoil  slopes  presents  some  problem  since  most  airfoil 
shapes  are  not  defined  by  simple  analytic  functions  but 
by  tabulated  points.  The  recommended  procedure  for  calcu- 
lating the  needed  slopes  is  to  perform  a cubic  spline  fit 
to  the  data.  However,  good  results  have  been  obtained  with 
STRANS  using  simpler  methods  such  as  linear  interpolation  to 
calculate  slopes.  The  effect  of  angle  of  attack  (ALPHA)  and 
flap  angle  (ALPHAF)  on  the  airfoil  shape  and  slope  functions 
are  treated  automatically  in  STRANS  using  the  specified 
input  quantities. 

UTRANS  automatically  treats  two  rigid  body  modes  of 
oscillation:  pitch  and  control  surface  rotation.  No  pro- 

gram modifications  are  required  as  the  different  airfoil 
boundary  conditions  are  determined  internally  using  the 
specified  input.  The  flag  for  one  mode  or  the  other  is 
the  hinge  point  (XH) ; if  XH=0  the  pitch  mode  is  treated  and 
if  XH/0  (0  _<  XH  _<  1)  control  surface  oscillations  are  con- 
sidered. The  method  is  not  restricted  to  these  modes  but 
would  require  modifications  to  FPU (I)  and  FPL (I)  in  subrou- 
tine INITAL  to  treat  various  flexible  modes  of  oscillation. 
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Grid  Design:  The  variable  mesh  size  technique,  used 

in  STRANS  and  UTRANS,  allows  considerable  flexibility  in 
the  choice  of  a computational  grid.  The  basic  reason  for 
this  capability  is  that  grid  points  may  be  concentrated  near 
the  airfoil,  where  a high  degree  of  resolution  is  desired. 

The  grid  spacing  may  then  be  expanded  away  from  the  air- 
foil out  to  the  boundaries  of  the  mesh,  thereby  providing 
for  an  economic  use  of  grid  points.  Economy  of  use  is  par- 
ticularly important  for  UTRANS  calculations  since  three 
values  (<J>°,  Re<J>:,  Iimf)1)  must  be  stored  at  each  nodal  point. 

The  basic  consideration  in  grid  selection  is  that 
grid  spacing  vary  smoothly  so  that  sudden  jumps  in  Ax  or 
Ay  should  be  avoided.  As  with  most  numerical  techniques, 
considerations  of  error  in  the  finite  difference  approx- 
imation indicate  that  the  ratio  of  the  spacing  between  ad- 
jacent grid  points  lie  between  1/2  and  2.  Some  additional 
guidelines  for  grid  design  are  that  grid  points  be  concen- 
trated near  the  leading  edge  to  resolve  the  blunt  stagna- 
tion point  singularity  and  the  pressure  singularity  in  a 
calculation  with  lift.  If  shocks  are  expected  in  the  cal- 
culation, grid  points  can  be  concentrated  near  their  expected 
location.  Choice  of  a grid  for  unsteady  calculations  fol- 
lows similar  guidelines  with  regard  to  concentration  of  grid 
points.  Unsteady  calculations  in  general  involve  lift  per- 
turbations to  the  steady  flow  so  that  the  leading  edge  or 
hinge  point  singularity  must  be  resolved.  Also,  if  the 
steady  flow  has  shocks  or  rapid  compressions,  wave  energy 
will  be  concentrated  there  so  that  a concentration  of  grid 
points  is  required  to  resolve  the  resulting  unsteady  pres- 
sure peak.  It  should  be  noted  that  UTRANS  requires  values 
of  0°,  generated  by  STRANS,  at  each  point  in  the  grid.  No 
provision  is  made  for  interpolating,  so  that  the  UTRANS  cal- 
culation must  have  the  same  basic  grid  as  the  STRANS  calcu- 
lation being  perturbed.  A provision  is  made  which  allows 
the  <£°  data  to  be  from  a more  refined  grid  which  was  obtained 
using  the  grid  halving  technique  on  the  equivalent  UTRANS 
grid.  The  utility  of  grid  halving  is  described  below. 

A basic  grid  which  has  provided  accurate  results 
in  calculations  performed  to  date  consists  of  about  50  grid 
lines  in  both  normal  and  streamwise  directions.  Approxi- 
mately half  of  the  grid  columns  (X  = constant)  are  concen- 
trated on  the  airfoil  and  the  grid  rows  (Y  = constant)  are 
divided  equally  above  and  below  the  airfoil.  Slight  vari- 
ations on  this  grid  are  given  in  the  test  cases  below.  Nu- 
merical studies  which  further  refined  the  basic  grid  des- 
cribed here,  have  shown  little  change  in  the  results  for 
pressure  coefficient  on  the  airfoils,  for  example,  with  the 
singular  exception  of  increased  resolution  of  a shock  wave. 
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It  is  judged  that  the  basic  grid  is  sufficiently  accurate 
for  most  engineering  applications. 

Farfield  Location  and  Update:  The  only  other  ques- 

tion in  grid  design  is  how  far  to  expand  the  grid  out  to 
the  farfield.  The  fact  that  the  perturbation  flowfield 
extends  to  large  lateral  distances  for  transonic  flows, 
makes  this  an  especially  important  question.  The  lateral 
extent  of  the  nonlinear  flowfield,  varies  considerably 
depending  upon  the  transonic  similarity  parameter,  so  that 
no  specific  answer  to  the  question  is  possible.  In  general 
though,  the  farfield  must  be  far  enough  away  from  the  air- 
foil that:  (i)  for  subsonic  freestream  flows,  the  asymp- 

totic expressions  used  for  setting  <j>ff  (Equation  37)  or 
<J>£f  (Equation  43)  be  an  accurate  representation  of  the  solu- 
tion, and  (ii)  for  supersonic  freestreams,  the  perturbation 
is  weak  and  the  flow  is  entirely  supersonic  on  the  boundary. 
For  subsonic  flow,  considerations  of  accuracy  of  the  analyt- 
ical farfield  expressions  indicate  that  Y ( JM) , Y(l)  = ±6  and 
X ( IM)  , x (1 ) = ±3  are  sufficient  at  least  for  <?  .95,  6 £ 
0.1,  a 2°  and  ft  and  K ^ 1.  It  is  noted  in  passing,  that 
numerical  studies  have  shown  that  the  farfield  approximation 
has  a weak  effect  on  steady  airfoil  pressures  and  a slightly 
stronger  effect  on  the  unsteady  pressures  on  the  airfoil. 

It  is  recalled  that  both  the  steady  and  unsteady  sub- 
sonic farfields  depend  upon  the  solution  through  the  nonlin- 
ear doublet  strength  (D0UBLT)  and  the  airfoil  circulation 
(GAMFF ) . These  quantities  are  thus  part  of  the  solution  pro- 
cess (as  described  in  Section  3.0)  and  couple  back  into  the 
solution  through  their  effect  on  the  farfield  and  wake.  The 
user  can  control  the  iteration  process  for  D0UBLT  and  GAMFF 
by  specifying  NGFF,  the  number  of  iterations  between  far- 
field updates  and  PGFF,  the  relaxation  parameter  for  updat- 
ing GAMFF.  It  is  suggested  that  NGFF  be  somewhat  greater 
than  the  number  of  grid  points  between  the  airfoil  and  the 
farfield  to  allow  the  updated  farfield  enough  iterations  to 
influence  the  entire  grid  before  updating  again.  A sugges- 
ted value  for  PGFF  is  1.5  which  has  the  effect  of  an  over- 
relaxation on  past  values  of  airfoil  circulation.  Such 
values  for  both  parameters  have  worked  well  in  the  past  al- 
though no  detailed  attempt  at  determining  optimum  values 
has  yet  been  made. 

Relaxation  Factors:  The  matter  of  the  choice  of 

optimum  relaxation  parameters  (cj)  is  a complex  one  and  has 
only  been  studied  in  detail  for  simple  cases.  In  the  case 
of  the  steady  transonic  potential  equation,  values  varying 
from  o)e  = 1.4  1.9  for  elliptic  points  and  = 0.75  -+  0.9 
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for  hyperbolic  points  have  been  used.  The  "best"  value 
can  vary  significantly  with  respect  to  the  transonic  sim- 
ilarity parameter  K,  however  limited  numerical  studies3 
indicate  a value  for  o)e  ^ 1.7  is  close  to  optimum  for 
slightly  supercritical  flows.  Similar  studies3  have  indi- 
cated that  the  relaxation  procedure  for  1 (UTRANS)  may 
become  unstable  for  values  of  much  greater  than  0.75 
and  the  best  values  for  and  seem  to  be  strongly  re- 
lated to  the  unsteady  similarity  parameters  K and  ft. 

Accuracy  and  Convergence:  Questions  of  accuracy 

and  convergence  of  the  numerical  schemes  implemented  in 
STRANS  and  UTRANS  are  of  considerable  importance  relative 
to  the  utility  of  the  programs  as  engineering  tools.  The 
accuracy  of  the  small  perturbation  approximation  for  steady 
flows  has  been  discussed  in  many  studies ( 4 > 5 > 6 ) and  STRANS 
results  are  consistent  with  these  findings3.  The  question 
of  the  accuracy  of  the  linearized  unsteady  perturbation 
approximation  used  in  UTRANS  is  an  open  one  but  results  to 
date3  have  been  quite  promising.  Of  particular  interest  to 
the  user  of  the  programs  is  the  degree  of  convergence  and 
grid  definition  required  for  practical  accuracy.  Clearly 
the  trade-off,  which  plagues  any  numerical  technique,  is 
the  cost  of  computer  time.  Some  comments  concerning  this 
question  are  therefore  appropriate. 

The  use  of  successive  grid  refinements  to  speed 
convergence  is  implemented  in  both  STRANS  and  UTRANS  in 
subroutine  D0UBLE  by  a simple  halving  of  the  grid  spacing. 
The  technique  has  been  found  to  be  quite  beneficial  to  the 
overall  efficiency  of  a calculation.  It  is  implemented  by 
starting  with  a coarse  grid  and  specifying  KEPS  > 1 and 
EPSGRD ( I ) , 1=1,  KEPS.  The  iteration  procedure  is  con- 
tinued on  the  coarse  grid  until  the  change  in  <j> 0 or  | <}> 1 | 
over  one  iteration  is  everywhere  less  than  EPSGRD (1),  a 
suggested  value  for  which  is  10  3 . A greater  degree  of 
convergence  on  an  interim  coarse  grid  seems  to  be  wasteful. 
The  grid  lines  are  then  doubled  and  the  iteration  proce- 
dure continued  until  the  change  in  <f>  is  less  than  EPSGRD  (2). 
If  this  is  the  final  grid  refinement  desired,  a suggested 
value  for  EPSGRD (2)  is  10“4.  Convergence  to  this  accuracy 
for  UTRANS  runs  requires  about  three  or  four  minutes  on  a 
CDC  6600  including  time  expended  on  both  coarse  and  refined 
grids.  STRANS  runs  may  be  of  comparable  or  in  general 
shorter  duration  depending  upon  similarity  parameters.  It 
is  believed  that  convergence  to  higher  accuracy  is  not  prof- 
itable since  changes  in  airfoil  pressure  distribution  are 
undetectable  (less  than  one  percent)  for  convergence  between 
10“ 4 and  10-5  for  example. 
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An  important  aspect  of  the  numerical  solution  pro- 
cedure is  the  iteration  process  for  determining  the  air- 
foil circulation.  This  procedure  and  the  corresponding 
treatment  of  the  Kutta  condition  and  updating  of  the  far- 
field  have  been  discussed  above.  Updating  the  circulation 
and  jump  in  the  potential  along  the  wake  has  a highly  non- 
linear effect  on  the  convergence  rate.  Typically,  (A<j>)m^x 
increases  by  an  order  of  magnitude  when  the  circulation  is 
updated.  Convergence  is  then  quite  rapid  until  the  next 
update  of  circulation  and  farfield  values.  Thus  an  addi- 
tional consideration  in  terminating  the  iteration  process, 
is  the  variation  in  circulation  from  one  update  to  the 
next.  The  change  in  circulation  (GAMFF)  over  NGFF  itera- 
tions as  well  as  the  change  in  the  nonlinear  doublet 
strength  (D0UBLT)  in  STRANS  calculations  are  good  indica- 
tors of  the  global  convergence  of  the  solution.  In  general, 
convergence  to  EPSGRD  of  10-4  is  equivalent  to  a change  in 
GAMFF  or  D0UBLT  over  NGFF  iterations  of  less  than  10“  . 

If  this  is  not  true,  it  is  suggested  that  the  calculation 
be  continued  to  a greater  degree  of  grid  convergence.  This 
will  insure  convergence  of  integral  quantities  such  as  air- 
foil force  coefficients  to  acceptable  engineering  accuracy. 

It  is  reiterated  that  the  comments  given  in  this  sec- 
tion are  intended  only  to  summarize  past  experience  with 
the  programs  and  the  suggested  values  for  various  parameters 
are  not  to  be  taken  as  hard  and  fast  rules.  Experimentation 
with  such  parameters  are  encouraged  and  in  fact  essential 
to  the  efficient  operation  of  the  programs. 


6 . 2 Summary  of  Recommended  Computational  Procedures 


Certain  standard  practices  used  in  any  large  scale 
numerical  computations  should  be  followed  in  exercising 
the  STRANS  and  UTRANS  programs.  For  instance,  when  initi- 
ating a series  of  calculations  for  a new  airfoil  it  is  ad- 
visable to  make  a short  run  of  perhaps  25  cycles.  In  this 
way  all  input  quantities  can  be  verified  and  possible  ad- 
justments to  the  grid  may  be  desired  based  on  the  progress 
of  the  calculation.  Such  a calculation  would  also  be  use- 
ful in  estimating  convergence  and  therefore  computer  time 
requirements  for  later  calculations.  The  peculiarities  of 
transonic  flow  computations  make  certain  additional  proced- 
ures advisable.  One  such  general  procedure  which  has  worked 
well  is  described  in  this  section.  It  makes  good  use  of 
the  restart  option  as  well  as  the  successive  grid  refine- 
ment technique. 
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It  is  noted  that  the  relaxation  procedure  for,  in 
this  case,  the  transonic  potential  equation  is  a systematic 
method  for  getting  from  one  solution  to  another  or  from  an 
initial  guess  to  the  desired  solution.  The  procedure  sum- 
marized here  thus  makes  use  of  this  inherent  characteris- 
tic to  generate  a matrix  of  approximate  solutions  for  var- 
ious Mach  numbers  and  angles  of  attack  for  a given  airfoil, 
which  can  be  used  as  relatively  good  initial  guesses  for 
more  refined  or  highly  convergent  solutions.  For  steady 
solutions  the  procedure  involves  initiating  a calculation 
for  a given  airfoil  on  a relatively  coarse  grid  (10  x 10 
to  25  x 25)  for  a subcritical  Mach  number  for  solutions 
with  Mo3  < 1 or  a completely  supersonic  Mach  number  if  solu- 

tions with  M^  > 1 are  desired.  In  either  case,  STRANS  uses 
the  linearized  subsonic  or  supersonic  solutions  respectively 
to  generate  an  initial  guess.  This  initial  case  is  relaxed 
to  a relatively  loose  degree  of  convergence  (VL0-3)  and  the 
restart  option  is  used  in  subsequent  runs  to  bootstrap  up 
or  down  in  Mach  number  if  subsonic  or  supersonic  solutions 
respectively  are  desired.  These  interim  solutions  can  then 
be  used  to  bootstrap  up  in  angle  of  attack.  Increments  in 
Mgch  number  of  about  .05  and  in  angle  of  attack  of  about 
1 have  generally  been  used  to  date  in  generating  these 
initial  solutions.  Any  one  of  these  calculations  can  later 
be  further  refined  using  the  grid  halving  technique  and 
ultimately  continued  to  the  recommended  10-1*  convergence 
level . 


The  procedure  for  unsteady  calculations  with  UTRANS 
is  quite  similar.  In  this  case,  the  procedure  is  initiated 
on  a coarse  grid  for  the  quasi-steady  or  k = 0 case,  using 
any  fully  refined  and  converged  steady  solution.  The  re- 
duced frequency  is  incrementally  increased  in  steps  of 
about  0.1  in  subsequent  calculations,  until  the  desired 
value  is  reached.  As  before  the  solution  can  then  be  re- 
fined and  continued  to  the  recommended  convergence  level. 

Clearly  many  variations  on  the  above  procedure  are 
possible  and  the  procedure  is  by  no  means  optimum  for  any 
given  calculation.  It  is  for  this  reason  that  the  syste- 
matic procedure  outlined  above  has  not  been  programmed  into 
STRANS  or  UTRANS. 
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7.0 


SAMPLE  CASES 


Detailed  input  and  sample  output  for  sequences  of 
STRANS  and  UTRANS  runs  are  presented  in  this  section.  The 
sample  cases  attempt  to  implement,  at  least  in  part,  the 
suggestions  presented  in  the  previous  section. 


7.1  STRANS  Test  Cases 


A sequence  of  computer  runs  are  described  in  this 
section,  which  calculate  the  steady  transonic  flow  over  a 
10  percent  thick,  symmetric  circular  arc  airfoil  for  the 
following  conditions: 


= 0.7 

00 


M = 0.8 

oo 


a = 0° 
a = 2° 


The  individual  runs  required  to  complete  these  cases  (i.e., 
achieve  convergence  oLO- 4 on  a refined  grid)  are  described 
in  the  run  log  given  in  Table  1.  The  table  lists  the  re- 
start tape  read  by  each  run,  the  grid(s)  used,  total  grid 
iterations,  convergence  achieved  on  each  grid  and  finally 
the  tape  dump  generated.  The  "coarse"  29x29  grid  on  which 
the  calculations  were  initiated  is  given  in  the  input  to 
Run  IS  below.  The  "refined"  54x55"  grid  is  generated  by 
halving  the  grid  spacing  as  described  earlier.  Runs  2S, 

5S,  6S,  7S  and  10S,  marked  by  * in  the  table,  are  believed 
to  be  of  engineering  accuracy  with  respect  to  grid  refine- 
ment and  convergence. 

The  runs  described  here  implement  two  approaches  to 
the  calculation  of  lifting  cases,  using  interim  solutions, 
and  it  is  worth  commenting  briefly  on  the  relative  effi- 
ciency of  the  approaches.  Runs  1,  3,  4 and  5,  for  = 0.7, 
generate  a solution  with  a = 2 0 by  increasing  a in  incre- 
ments using  the  interim  coarse  grid  solutions.  That  is. 

Run  3 uses  file  IS. a as  its  initial  guess  for  a = 1°  and 
Run  4 uses  file  3S  as  its  initial  guess  for  a = 2°.  After 
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achieving  a certain  degree  of  convergence  on  the  coarse 
grid  Run  4 refines  the  grid  and  the  solution  is  completed 
in  Run  5.  The  sequence  of  = 0.8  runs  (Runs  7,  8,  9,  10) 
on  the  other  hand  use,  in  each  case,  the  refined  grid  solu- 
tion as  an  initial  guess  for  the  successively  larger  angles 
of  attack.  Accounting  for  the  fact  that  iterations  on  the 
refined  grid  require  about  a factor  of  4 more  computing 
time,  it  would  geem  that  the  sequence  of  runs  for  the 

= 0.7,  a = 2 solution  required  almost  a factor  of  4 less 
computing  time  than  the  = 0.8,  a = 2°  solution.  The  rea- 
son for  the  increased  efficiency  is  that  the  convergence  of 
lifting  cases  is  generally  controlled  by  the  iteration  pro- 
cess for  determining  the  circulation.  Thus  it  would  seem 
to  be  most  efficient  to  achieve  relatively  good  convergence 
for  the  circulation  before  proceeding  to  the  refined  grid. 
This  of  course  neglects  the  fact  that  supercritical  cases 
such  as  = 0.8  generally  require  greater  computing  effort 
than  subcritical  (M^  = 0.7)  cases.  However  the  increased 
effort  for  supercritical  cases  is  not  a factor  of  4 by  any 
means  so  that  this  comparison  would  seem  to  justify  the  con- 
clusion given.  The  comparison  of  "bootstrapping"  techniques 
given  here,  also  reinforces  the  point  made  earlier  that  the 
user  can  effect  considerable  savings  by  exercising  some 
judgment  in  the  way  calculations  are  performed. 


Table  1.  Sequence  of  Runs  for  STRANS  Sample  Cases  (*  Converged  Solution) 


Run 

M*> 

a 

Restart 

Tape 

Used 

Grid 

Type 

Grid 

Iterations 

Convergence 

Achieved 

Tape  Dump 
Generated 

IS 

.7 

0° 

— 

29*29 

54*55 

58 

ro  ro 
1 1 
o o 

H H 

IS.  a 
lS.b 

* 2S 

.7 

0° 

lS.b 

54*55 

85 

10'4 

2S 

3S 

.7 

1° 

IS. a 

29*29 

44 

10-3 

3S 

4S 

.7 

2° 

3S 

29*29 

54*55 

92 

10-3 
10  J 

4S.a 

4S.b 

* 5S 

.7 

2° 

4S.b 

54*55 

100 

1.2*10"4 

5S 

* 6S 

.75 

0° 

2S 

54*55 

68 

10-4 

6S 

* 7S 

.8 

0° 

6S 

54*55 

75 

10-4 

7S 

8S 

.8 

1° 

7S 

54*55 

100 

10-3 

8S 

9S 

.8 

2° 

8S 

54*55 

100 

1 . 6*10-3 

9S 

*10S 

.8 

2° 

9S 

54*55 

200 

4.7*10'4 

10S 
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7.1.1 


Input  for  STRANS  Test  Cases 


The  card  input  for  each  of  the  STRANS  runs  described 
above  is  given  in  this  section: 


• Run  IS:  no  tape  read,  generate  files  IS. a,  lS.b. 


***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=0 , 

$END 

$IN 

X ( 1 ) =-4 . , -3. 2, -2. 4, -1. 6,-1. ,-.6, -.35, -.2, -.075,0. , .05,  .11, 
.18,. 26,. 36,. 5,. 64,. 74,. 82,. 88., 94,1., 1.1, 1.25, 1.45, 1.75, 
2.25,3. ,3.75, 

Y(l)=-7. , -5. 4, -3. 8, -2. 7, -1.9, -1.4, -1.1, -.85, -.65, -.5, -.36, 
-.23, -.1366, -.05,0.,. 05,. 13 66,. 23,. 36,. 5,. 65,. 85, 1.1, 1.4, 
1. 9, 2. 7, 3. 8, 5. 4, 7.  , 

IM=2  9 , 

JM=  2 9 , 

ILE=10 , 

ITE=22 , 

JW= 1 5 , 

M8= . 7 , 

GAM= 1.4, 

DEL=. 1, 

ALPHA=0 . , 

GAMFF=0 . , 

0MEGAH= .75, 

0MEGAE=1 . 5 , 

0MEGAP= .75, 

EPSC0L=5 . E-5 , 

EPSGRD (1) =1 . E-3 , 1 . E-3 , 

XH=0 . , 

NDUMP=2000 , 

NC0L=1O , 

NGRID=150 , 

NGFF=10 , 

PGFF=1 . 5 , 

KEPS-2  , 

NPRINT=10, 

ALPHAF=0 . , 

$END 
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Run  2S:  read  file  lS.b,  generate  file  2S. 


***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=1 , 

$END 

$IN 

0MEGAH=.75, 
0MEGAE=1 . 7 , 
0MEGAP=  .75, 
EPSC0L=5 . E-5 , 
EPSGRD (1 ) =1 . E-4 , 
NDUMP=2000, 
NC0L=1O, 
NGRID=100 , 
NGFF=10 , 

PGFF=1. 5, 

KEPS=1 , 

NPRINT=10 , 

$END 


• Run  3S : read  file  IS. a,  generate  file  3S. 


***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=1 , 

$END 

$IN 

ALPHA=. 0174532925, 
GAMFF=1 . , 

0MEGAH= .75, 
0MEGAE=1 . 7 , 
0MEGAP=. 75, 
EPSC01i=5 . E-5 , 
EPSGRD (l)=l.E-3, 
NDUMP=2000, 
NC0(L=1O, 

NGRID=100 , 

NGFF=1 0 , 

PGFF=1. 5, 

KEPS=1, 

NPRINT=10 , 

IK=1, 

$END 


-55- 


Run  4S:  read  file  3S,  generate  files  4S.a  and  4S.b 


***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=1 , 

$END 

$IN 

ALPHA= . 034906585, 
GAMFF=1. , 

0MEGAH=. 75, 

0MEGAE= 1.7, 

0MEGAP= -75, 

EPSC0L=5 . E-5 , 

EPSGRD ( 1 ) =1 . E- 3 , 1 . E- 3 , 
NDUMP=2  000 , 

NC0L=1O , 

NGRID=150 , 

NGFF=10 , 

PGFF=1. 5, 

KEPS=2 , 

NPRINT=10 , 

IK=1 , 

$END 


• Run  5S : read  file  4S.b,  generate  file  5S. 


***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=1 , 

$END 

$IN 

0MECAH=  .75, 
0MEGAE=1. 7, 
0EMGAP= .75, 
EPSC0L=5 . E-5 , 
EPSGRD (l)=l.E-4, 
NDUMP=2  000 , 
NC0L=1O , 

NGRID=1 00 , 
NGFF=10 , 

PGFF=1 .5, 

KEPS=1 , 

NPRINT=10 , 

$END 
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Run  6S:  read  file  2S,  generate  file  6S. 


***  CIRCULAR  ARC  *** 


$C0NTRL 
I TAPE =1 , 

$END 

$IN 

M8=. 75, 

0MEGAH= .75, 
0MEGAE=1 . 7 , 
0MEGAP=. 75, 
EPSC0L=  5 . E- 5 , 
EPSGRD ( 1 ) =1 . E-4 , 
NDUMP=2000, 
NC0L=1O , 
NGRID==100 , 
NGFF=10 , 

PGFF=1. 5, 

KEPS=1 , 
NPRINT=10, 

IK=1 , 

$END 


• Run  7 S : read  file  6S,  generate  file  7S. 


***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=1 , 

$END 

$IN 

M8=. 8, 

0MEGAH= . 75, 
0MEGAE=1 . 7 , 
0MEGAP= . 75 , 
EPSC0L=  5 . E- 5 , 
EPSGRD (1)=1. E-4, 
NDUMP=2000, 
NC0L=1O , 
NGRID=100 , 
NGFF=10 , 

PGFF=1. 5, 

KEPS=1 , 

NPRINT:=10 , 

IK-1, 

$END 
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Run  8S:  read  file  7S,  generate  file  8S 


***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=1 , 

$END 

$IN 

ALPHA=. 0174532925, 
GAMFF=1 . , 

0MEGAH= .75, 
0MEGAE=1 . 7 , 
0MEGAP=. 75, 
EPSC0L=5 . E-5 , 
EPSGRD (1) =1 .E-4 , 
NDUMP=2000, 
NC0L=1O, 

NGRID=10  0 , 

NGFF=10 , 

PGFF=1 . 5 , 

KEPS=1 , 

NPRINT=10 , 

IK=1 , 

$END 


• Run  9S:  read  file  8S,  generate  file  9S 


***  CIRCULAR  ARC  *** 

$C0NTRL 
ITAPE=1 , 

$END 

$IN 

ALPHA= . 034906585, 

GAMFF=1.  , 

0MEGAH=  .75, 

0MEGAE=1.7, 

0MEGAP= .75, 

EPSC<?L=5.E-5, 

EPSGRD ( 1 ) =1 . E-4 , 

NDUMP=2000 , 

NC<2tL=10, 

NGRID=10  0 , 

NGFF=10 , 

PGFF=1. 5, 

KEPS=1 , 

NPRINT=10, 

IK=1, 

$END 
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Run  10S:  read  file  9S,  generate  file  10S. 


***  CIRCULAR  ARC  *** 


$C0NTRL 
I TAPE =1 , 

$END 

$IN 

0MEGAH= .75, 
0MEGAE=1.7, 
0MEGAP=.75, 
EPSC?L=5.E-5, 
EPSGRD ( 1 ) =1 . E-4  , 
NDUMP=2  000 , 
NC0L=1O, 
NGRID=2000 , 
NGFF=10 , 

PGFF=1. 5, 

KEPS=1, 

NPRINT=10, 

$END 


7.1.2  Sample  Output  for  STRANS  Test  Cases 


The  following  pages  contain  a 
uous  commentary  output  for  the  first 
addition  to  the  final  printed  output 


sample  of  the  contin- 
9 cycles  of  run  IS  in 
for  all  runs. 
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TERATION  1 COLUMN  9 FAILED  TO  CONVERGE  ERR  4 .40315E-03  Jf  11 


• Sample  Output  from  Run  IS 
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HINGE  POINT  = 0. 
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CIRCULAR  ARC 
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THICKNESS  PATIO  = .10000E*00 

AIRFOIL  ANGLE  OF  ATTACK  (RAOIANS)  = 0. 
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7.2 


UTRANS  Test  Cases 


A sequence  of  UTRANS  runs  are  described  in  this  sec- 
tion, which  calculate  the  unsteady  flow  perturbation  for 
the  circular  arc  airfoil  oscillating  in  angle  of  attack 
at  the  following  conditions: 


M = 0.7  , a = 

oo  ' 


Mm  = 0.8  , a = 

OO 


The  individual  runs  required  to  complete  these  cases  (i.e., 
achieve  convergence  of  ^10" 4 on  a refined  grid)  are  des- 
cribed in  the  run  log  in  Table  2.  The  table  presents  the 
mean  flow  conditions  (M^, a0 ) and  reduced  frequency  (k) , the 
tape  dump  of  the  steady  solution  and  restart  tape  read  for 
each  run,  the  grid(s)  used,  total  grid  iterations,  conver- 
gence achieved  on  each  grid  and  finally  the  tape  dump  gen- 
erated. All  runs  are  on  the  same  coarse  and  refined  grids 
described  above.  Runs  2U,  5U,  7U  and  10U,  marked  by  * in 
the  table,  are  the  final  converged  results  for  the  cases 
listed  above.  The  runs  were  performed  in  the  order  shown 
and  implement  the  "bootstrap"  technique  for  getting  from 
one  unsteady  solution  to  the  other  as  described  in  the 
previous  section.  The  required  input  for  each  run  and 
sample  output  are  now  presented. 


( k 

>° 

( k 

,o! 


k = 


k = 


0.0 

0.2 

0.0 

0.2 
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Table  2.  Sequence  of  Runs  for  UTRANS  Sample  Cases  (*  Converged  Solution) 


Run 

M» 

a 

k 

STRANS 

Tape 

Restart 

Tape 

Grid 

Type 

Grid 

Iterations 

Convergence 

Achieved 

Tape  Dump 
Generated 

1U 

.7 

0° 

0.0 

2S 

— 

29x29 

54x55 

45 

3xl0-3 

10"J 

lU.a 

lU.b 

* 2U 

.7 

0° 

0.0 

2S 

lU.b 

54x55 

200 

1.85xl0-4 

2U 

3U 

.7 

0° 

0.1 

2S 

lU.a 

29x29 

127 

10“3 

3U 

4U 

.7 

0° 

0.2 

2S 

3U 

29x29 

54x55 

136 

10-3 
10  J 

4U.a 

4U.b 

D 

i 

.7 

0° 

0.2 

2S 

4U.b 

54x55 

200 

1C"4 

5U 

6U 

.8 

0° 

0.0 

7S 

lU.a 

29x29 

54x55 

189 

10-3 
10  J 

6U.a 

6U.b 

* 7U 

.8 

0° 

0.0 

IS 

6U.b 

54x55 

199 

IQ'4 

7U 

8U 

.8 

0° 

0.1 

IS 

6U.a 

29x29 

187 

10~3 

8U 

9U 

.8 

0° 

0.2 

IS 

8U 

29x29 

54x55 

175 

10-3 

10-3 

9U.a 

9U.b 

*10U 

.8 

0° 

0.2 

IS 

9U.b 

54x55 

200 

ID"4 

10U 
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7.2.1 


Input  for  UTRANS  Test  Cases 


The  card  input  for  each  of  the  UTRANS  runs  described 
above  is  given  in  this  section. 


• Run  1U:  read  file  2S,  no  restart  tape  read;  gen- 

erate files  lU.a  and  lU.b. 

***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=0 , 

$END 
$ IN 

X(l)=-4. ,-3.2, -2. 4,-1. 6,-1. , -.6, -.35,-. 2, -.075,0. ,.05,. 11, 
.18, .26, .36, .5, .64, .74, .82, .88, .94,1. ,1.1,1.25,1.45,1.75, 
2.25,3. ,3.75, 

Y(l)=-7. ,-5.4, -3. 8, -2. 7,-1. 9, -1.4, -1.1, -.8 5, -.65, -.5, -.36, 
-.23, -.1366, -.05,0. , .05, .1366, .23,. 36,. 5,. 65,. 85, 1.1, 1.4, 
1.9, 2. 7, 3. 8, 5. 4, 7.  , 

GAMFF= ( 2 . , 0 . ) , 

IM=29 , 

JM=29 , 

ILE=10 , 

ITE=22 , 

JW=1 5 , 

0MEGAH= .75, 

0MEGAE=1. 7, 

0MEGAP= . 75 , 

EPSGRD ( 1 ) =3 . E- 3 , 1 . E-3 , 

NDUMP=2000 , 

NGRID=200 , 

NGFF=10 , 

PGFF=1 . 5 , 

KEPS=2 , 

NPRINT=10 , 

SMALLK=0. , 

XH=0 . , 

$END 
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• Run  2U:  read  file  2S,  restart  file  lU.b;  generate 

file  2U. 


***  CIRCULAR  ARC  *** 


$C0NTRL 

ITAPE=1, 

$END 

$IN 

0MEGAH= .75, 
0MEGAE=1 . 7 , 
0MEGAP= .75, 
EPSGRD ( 1 ) =1 . E-4 , 
NDUMP=2000 , 
NGRID=2  0 0 , 
NGFF=10 , 

PGFF=1. 5, 

KEPS=1 , 

NPRINT=10 , 

$END 


• Run  3U : read  file  2S,  restart  file  lU.a;  generate 

file  3U. 


***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=1 , 

$END 

$IN 

GAMFF= (2. 3,-1. ) , 
SMALLK=. 1, 
0MEGAH= .75, 
0MEGAE=1 . 7 , 
0MEGAP=. 75, 
EPSGRD ( 1 ) =1 . E- 3 , 
NDUMP=2  000 , 
NGRID=200, 
NGFF=10 , 

PGFF=1. 5, 

KEPS=1 , 

NPRINT=10 , 

IK=1 , 

$ END 
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Run  4U : read  file  2S,  restart  file  3U,  generate 

file  4U . 


***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=1 , 

$END 
$ IN 

GAMFF= (2 . ,-.5) , 

SMALLK= . 2 , 

0MEGAH= .75, 

0MEGAE=1 . 7 , 

0MEGAP= .75, 

EPSGRD ( 1 ) =1 . E-3 , 1 . E-3 , 
NDUMP=2000 , 

NGRID=200 , 

NGFF=10 , 

PGFF=1 . 5 , 

KEPS=2 , 

NPRINT=10, 

IK=1 , 

$END 


• Run  5U:  read  file  2S,  restart  file  4U;  generate 

file  5U. 


***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE==1 , 

$END 
$ IN 

0MEGAH= .75, 
0MEGAE=1 . 7 , 
0MEGAP= .75, 
EPSGRD (l)=l.E-4, 
NDUMP=2000 , 
NGRID=2  00 , 
NGFF=10 , 

PGFF=1 . 5 , 

KEPS=1 , 

NPRINT=10 , 

$END 
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• Run  6U : read  file  IS,  restart  file  lU.a;  generate 

files  6U.a,  6U.b. 


***  CIRCULAR  ARC  *** 

$C0NTRL 

ITAPE=1, 

$END 

$IN 

GAMFF= (2.6,0.)  , 

0MEGAH= .75, 

0MEGAE=1. 7, 

0MEGAP= .75, 

EPSGRD ( 1 ) =1 . E-3 , 1 . E-3 , 

NDUMP=2000 , 

NGRID=200 , 

NGFF=10 , 

PGFF=1. 5, 

KEPS=2 , 

NPRINT=10, 

IK=1 , 

$END 


• Run  7U:  read  file  IS,  restart  file  6U.b;  generate 

file  7U. 


***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=1 , 

$END 

$IN 

0MEGAH= .75, 
0MEGAE=1 . 7 , 
0MEGAP= .75, 
EPSGRD (l)=l.E-4, 
NDUMP=2000 , 
NGRID=200 , 
NGFF=10 , 

PGFF=1 . 5 , 

KEPS=1 , 

NPRINT=10 , 

$END 
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Run  8U:  read  file  7S,  restart  file  6U.a;  generate 

file  8U. 


***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=1 , 

$ END 
$ IN 

GAMFF= ( 3 . 8 , -1 . ) , 
SMALLK=. 1, 
0MEGAH=. 75, 
0MEGAE=1 . 7 , 
0MEGAP= .75, 
EPSGRD ( 1 ) =1 . E-3 , 
NDUMP=2000 , 
NGRID=100 , 
NGFF=10 , 

PGFF= 1 . 5 , 

KEPS=1 , 

NPRINT=10 , 

IK=1 , 

$END 


• Run  9U:  read  file  7S,  restart  file  8U;  generate 

files  9U.a  and  9U.b. 

***  CIRCULAR  ARC  *** 


$C0NTRL 
ITAPE=1 , 

$END 

$IN 

GAMFF= (3. 5,-1. 5) , 
SMALLK=.  2, 

0MEGAH=. 75, 

0MEGAE=1 . 7 , 

0MEGAP=. 75, 

EPSGRD (1 )=1. E-3, 1. E-3, 
NDUMP=2000 , 

NGRID=2  00 , 

NGFF=1 0 , 

PGFF=1. 5, 

KEPS=2 , 

NPRINT=10 , 

IK=1 , 

$END 
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# Run  10U:  read  file  7S,  restart  file  9U.b;  generate 

file  10U. 

***  CIRCULAR  ARC  *** 


$ C0NTRL 
ITAPE=1 , 

$END 

$IN 

0MEGAH= .75, 
0MEGAE=1. 7, 
0MEGAP=. 75, 
EPSGRD ( 1 ) = 1 . E- 4 , 
NDUMP=2000 , 
NGRID=2  0 0, 
NGFF=10 , 

PGFF=1. 5, 

KEPS= 1 , 

NPRINT=10 , 

$END 


7.2.2  Sample  Output  for  UTRANS  Test  Cases 


The  following  pages  contain  a sample  of  the  contin- 
uous commentary  output  for  the  first  19  cycles  of  run  1U 
in  addition  to  the  final  printed  output  for  all  runs. 
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• Sample  Output,  for  Run  1U 
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APPENDIX  A 


FORTRAN  LISTING  OF  STRANS 


A FORTRAN  listing  of  the  source  deck  for  the  STRANS 
program  is  presented  in  the  following  pages.  The  program, 
as  configured  here,  requires  638K  words  to  load  and  518K 
words  to  execute. 
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PROGRAM  STRANS  ( INPUT , OUTPUT  *TAPE5= INPUT » TAPE6=0UTPUT  * TAPE 7* 

1  TAPES) 

REAL  K »M8 

COMMON  /DELTA/  OX (99) ,DY (99) , AX1 (99) »AX2 (99) »BXl (99) ,BX2 (99) , 

1 CX (99) *AY1 (99) ,AY2(99>  * X ( 1 00 ) ♦ Y ( 1 00 ) ♦FPU ( 99) ♦ FPL (99) fPHIUB(99) ♦ 

2 IM,IM1,JM,JM1,JW»JWP1,JWM1,ITE,ILE»DYBU1,DYBU2,DYBL1,DYBL2, 

3 ALPHA ♦ DEL »M8»GAM»K  »NDB  »XH»TITLE(8) »DOUB»DOUBLT  »ALPHAF 
COMMON  /COE FF/  A (99) ,B (99) ,C (99) ,D (99) ,PHI ( 10000) 

COMMON  /GAMMA/  GAmTEI »GAMTE»PGFF  »GAMFF 

DIMENSION  PHIOC (99) tPHlOG(99) ♦ OMEGA ( 99 ) ♦ V < 99 ) tEPSGRD (3) 

NAMELIST  /IN/  X ♦ Y tM8 ,GAM ♦ DEL ♦ ALPHA ♦GAMFF ♦ IM, JM ♦ ILE ♦ I TE ♦ JW ♦ 

1 OMEGAH  »OMEGAE  »OMEGAP ,EPSCOL ,EPSGRD ♦ XH  »NDIJMP  ,NCOL  tNGR ID ♦ 

2 NGFFtPGFFtKEPStlKtNPRINTtALPHAF 
NAMELIST  /CONTRL/  ITAPE 

C 

C TO  RESTART  PROGRAMt  COPY  THE  DATA  FOR  RESTART  FROM  T APE7  TO  A DISC 
C FILE  TAPE8,  POSITION  TAPE7  AT  THE  END  OF  THE  LAST  FILE  ON  THE  TAPE 
C SO  NEW  DATA  MAY  BE  WRITTEN  ON  THE  TAPE  WITHOUT  LOSING  ANY  OF  THE 
C OLD  DATA 
C 

READ  (5t912)  (TITLE ( I > • 1=1 »8) 

READ  ( 5 tCONTRL ) 

IF  (ITAPE. EQ.O)  GO  TO  10 
C READ  DATA  FROM  RESTART  TAPE 

READ  (8)  NITERG»IM»IM1,JM,JM1,JW,JWP1,JWM1,ITE»ILE»K,DEL, 

1 ALPHA, GAMTEI ,GAMFF »NDB , XH,M8 ,GAM , DYBU1 ,DYBU2 » DYBL 1 »DYBL2 » 

2 DOUBtDOUBLT  »ALPHAF 

READ  (8)  (X(I> ,I=1,IM) 

READ  (8)  (Y(I) ,1=1 ♦ JM) 

READ  (8)  (FPU (I) ,FPL(I) ,PHIUB(I) »I=ILE»ITE) 

READ  (8)  (AX1 (I) ,AX2(I) ,BX1 (I) ,BX2(I) ,CX(I) ,DX(I) ,I=1,IM1) 

READ  (8)  (AY1 (I) ,AY2(I) ,DY (I) ♦ 1=1, JM1) 

L=IM*JM 

READ  (8)  (PHI ( I) ♦ 1=1 ,L) 

I K = 0 

READ  (5, IN) 

C THE  IK  OPTION  IS  USED  TO  BOOT  STRAP  TO  DIFFERENT  MACH  NUMBERS, 

C AIRFOIL  THICKNESSES  AND/OR  ANGLES  OF  ATTACK 
C MAKE  SURE  YOU  HAVE  INPUT  THE  NEW  M8,  DEL  AND/OR  ALPHA 
IF  (IK. EQ.O)  GO  TO  A 

K= ( 1 . -M8**2) / ( ( 1 . ,GAM) *DEL*M8**2) «* .6666666667 
CALL  INITAL 
CALL  FARFLD 
4 CONTINUE 

SK=SQRT ( ABS (K ) ) 

GAMTE=GAMTE 1 
WRITE  (6,900) 

WRITE  (6,901)  NITERG 

NITERG=0 

GO  TO  15 

C START  PROBLEM  FROM  SCRATCH 
10  CONTINUE 
READ  (5, IN) 

K=  ( 1 .-M8**2)/ ( ( 1 .,GAM) *DEL*M8**2) **.6666666667 
SK=SQRT (ABS(K) ) 
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GAMTE=GAMEE 

GAHTE 1=GAMFF 

NITERG=0 

NDB=0 

IM1=IM-1 

JH1=JM-1 

JWP1=JW*1 

JWM 1 = JW- 1 

IMJM=IM*JM 

DO  5 L = 1 * I M JM 

PHI (L)=0. 

5 CONTINUE 

C INITIALIZE  FINITE  DIFFERENCE  COEFFICIENTS  AND  FAREIELO 
CALL  INITAL 
CALL  FARELO 
IF  (K.GT.O.)  GO  TO  45 

C INITIAL  GUESS  FOR  SUPERSONIC  CASE  (INTERIOR  ONLY) 

DO  40  1=3. IM1 
M = ( I — 1 ) •JM 
DO  41  J=2» JM1 
L = M>  J 
PHI ( L ) =0 • 

SKY=SK*ABS(Y(J> ) 

IF  < X < I ) .LT.SKY.OR.X(I) .GT.SKYM.)  GO  TO  41 
IE  (Y(J).GE.O.)  PHI (L)=-EU(X(I)-SKY)/SK 
IF  (Y(J).LT.O.)  PHI (L)=-EL(X(I)-SKY)/SK 
41  CONTINUE 
40  CONTINUE 
GO  TO  46 

45  CONTINUE 

C INITIAL  GUESS  FOR  SUBSONIC  CASE  (INTERIOR  ONLY) 
CON=DOUB/ (3. 14 159265* SK) 

DO  6 1=3. IM1 
M=(I-1)*JM 
DO  7 J=2*JM1 
L=M*  J 

IF  (X(I).EQ.O..AND.Y(J).EQ.O.)  GO  TO  8 
PHI (L)=C0N*X(I)/(X(I)**2*K*Y(J)**2> 

IF  (ABS(PHI (L) ) .GT.l.)  PHI (L ) =S IGN ( 1 . ♦ X ( I ) ) 

GO  TO  7 
8 CONTINUE 

PHI (L ) =PHI (L-JM) 

7 CONTINUE 

6 CONTINUE 

46  CONTINUE 
L=(ILE-2)*JM*JW 
PHIUB(ILE)=PHI (L) 

15  CONTINUE 

WRITE  (6* IN) 

WRITE  (6.900) 

CPCPB=DEL**. 6666666667/ ( ( 1 . ♦GAM) *M8**2) **.3333333333 
WRITE  (6.913)  K.CPCPB 
KGRD= 1 
50  CONTINUE 
ERROR=0. 
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N I T=N I TERG 
NITERG=NITERGM 

IF  (H0D(NITERG»NPRINT) .EQ.O)  CALL  PRINT (NIT) 

IF  (K.LT.O.)  GO  TO  51 
IF  (MOD(NITERGtNGFF) .NE.O)  GO  TO  51 
CALL  GAMFUN 
CALL  FARFLD 

WRITE  (6*910)  NITERG.GAMFF*GAMTE*DOUBLT 
51  CONTINUE 

DO  1 I*2*JM1 
V ( I ) =K 

1 CONTINUE 

C BEGIN  LOOP  ON  GRID 
DO  100  I = 3 » IM 1 
C CHECK  FOR  AIR  FOIL 
IFLAG=0 

IF  (ILE.LE.I.AND.I.LE.ITE)  IFLAG= 1 
M = ( 1-1 )*JM 

C SAVE  THIS  COLUMN  OF  PHI 
DO  2 J=2  * JM 1 
L=M*J 

PH  I OG ( J ) =PH I (L) 

2 CONTINUE 

C LOOP  BACK  POINT  FOR  COLUMN  ITERATION 

niterc=o 

150  CONTINUE 

NITERC=NITERC*1 
IF  (NITERC.GT.NCOL)  GO  TO  294 
C SAVE  PREVIOUS  PHI  FOR  COLUMN  ITERATION 
DO  3 J=2»JM1 
L = M*  J 

PHIOC ( J) =PHI (L) 

3 CONTINUE 

C BEGIN  LOOP  ON  COLUMN 
DO  160  J=2  * JM 1 
C CALCULATE  CELL  INDICES 
L=M  ♦ J 
LP=L* JM 
LL=L-JM 
LLL=LL-JM 
LB=L-1 
LT=L* 1 

TPHIR=PHI (LR) 

TPHIL=PHI (LL) 

IF  ( I .EO. ILE-1 .AND.J.EQ. JW)  PH  I (LR ) = .5* (PH  I (L* JM ) ♦PH IUB ( ILE ) ) 
IF  (I.EQ.ITEM.AND.J.EQ.JW)  PHI  (LL)  = .5*  (PHIUB  ( ITE)  *PHI  (L-JM)  ) 
C SET  UP  TRIDIAGONAL  MATRIX  TO  SOLVE  FOR  PHI(I*J) 

C A * PH  I ( I * J* 1 ) ♦ B * PHI ( I * J)  ♦ C * PH 1(1* J- 1 ) = D 
IF  (IFLAG.EQ.1.AND.J.EQ.JWP1)  GO  TO  250 
IF  (IFLAG.EQ.1.AND.J.EO.JWM1)  GO  TO  270 
IF  (IFLAG.EQ.l.AND.J.EQ.JW)  GO  TO  290 
PART=0 • 

IF  (I.LE.ITE)  GO  TO  201 
C KUTTA  CONDITION 

SIGI=(X(I)-1.)* (GAMFF-GAMTE ) /(X(IM1)-1.)*GAMTE 
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IF  (J.EQ.JWM1)  PART=.5*AY1 ( J)*SIGI 
IF  (J.EQ.JW)  PART=.5*(AY1 (J)-AY2(J> )*SIGI 
IF  (J.EQ.JWP1)  PART=-.5*AY2( J)*SIGI 

201  CONTINUE 

VV=K-AXl (I)* (PHI (LR)-PHI (L) ) -AX2 ( I ) * ( PHI (L ) -PHI (LL ) ) 

IF  (VV.LT.O.)  GO  TO  220 
C ELLIPTIC 

OMEGA (J)=OMEGAE 

B(J)=-(VV*(BX1 (I) ♦8X2(1) ) ♦ A Y 1 (J)*AY2(J> ) 

D ( J ) = -VV* ( BX 1 ( I ) *PH I (LR > *BX2 ( I ) *PH I (LL) ) »PART 
IF  (J.EQ.2)  GO  TO  202 
IF  (J.EO.JM1)  GO  TO  203 
At J)=AYI ( J) 

C ( J) =AY2 ( J) 

GO  TO  200 

C BOTTOM  BOUNDARY 

202  CONTINUE 

At J)=AY1 (J) 

D(J)=D(J)-AY2(J)*PHI (LB) 

GO  TO  200 
C TOP  BOUNDARY 

203  CONTINUE 

C ( J) =AY2 ( J) 

D(J)=D( J)-AY1 ( J)*PHI (LT) 

GO  TO  200 
220  CONTINUE 

IF  (V(J) .GT.O.)  GO  TO  240 
C HYPERBOLIC 

OMEGA (J) =OMEGAH 

VV=K-CX(I-1)* (PHI (L)-PHI (LL) ) -CX ( I -2 ) * ( PH  I (LL)-PHI (LLL) ) 
B(J)=VV*BX1 (I-l)-AYl ( J) -AY2 ( J) 

D(J)=VV*(BX1 ( 1-1 ) *PHI (LL)*BX2(I-1)*(PHI (LL)-PHI (LLL ) ) ) *PART 
IF  (J.EQ.2)  GO  TO  222 
IF  (J.EQ.JMi)  GO  TO  223 
A ( J) =AY1 ( J) 

C ( J ) =A Y2 ( J ) 

GO  TO  200 

C BOTTOM  BOUNDARY 

222  CONTINUE 
A(J)=AY1 (J) 

IF  (K.LT.O.)  GO  TO  224 
D { J) =D ( J) -AY2 ( J) *PHI (LB) 

GO  TO  200 

224  CONTINUE 

B(J)=B(J)-2.*SK*CX(I-1)/DY(2) 

D(J)=D<J)-2.*SK*(CX(I-1>*PHI (LL) -CX ( 1-2) * (PHI (LL)-PHI (LLL) ) )/DY(2) 
GO  TO  200 
C TOP  BOUNOARY 

223  CONTINUE 

C ( J) =AY2 ( J) 

IF  (K.LT.O.)  GO  TO  225 
D< J)=D(J)-AY1 ( J ) *PH I (LT) 

GO  TO  200 

225  CONTINUE 

B< J)=B(J)-2.*SK*CX (I -I )/DY( J-l ) 
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D(J)=D(J)-2.*SK*(CX(I-1)*PHI (LL)-CX(I-2)*(PHI (LL)-PHI (LLL ) ) )/ 
1 DY(J-l) 

GO  TO  200 
C PARABOLIC 
240  CONTINUE 

OMEGA (J)=OMEGAP 

B(J)=VV*BXl (I-l)-AYl (J)-AY2(J) 

D(J)=VV*(BX1  (1-1)48X2  (I- 1) >*PHI  (LL ) -VV«BX2 ( I- 1 > *PHI (LLL ) »PART 
IF  (J.E0.2)  GO  TO  242 
IF  (J.EO.JM1)  GO  TO  243 
A ( J) =AYI ( J) 

C ( J ) = A Y2 ( J ) 

GO  TO  200 

C BOTTOM  BOUNOARY 

242  CONTINUE 

A ( J) =AY1 (J) 

D(J)=D(J)-AY2(J)*PHI (LB) 

GO  TO  200 
C TOP  BOUNDARY 

243  CONTINUE 

C ( J) =AY2 ( J) 

D ( J) =D ( J) -AY1 ( J)*PHI (LT) 

GO  TO  200 

C BODY  BOUNDARY  ABOVE 
250  CONTINUE 

W = K-AX1  ( I ) * ( PHI  (LR)-PHI  (L)  ) -AX2  ( I ) * (PHI  (L)-PHI  (LL)  > 

IF  (VV.LT.O.)  GO  TO  260 
C ELLIPTIC 

OMEGA (J)=OMEGAE 
A ( J) =DYBU1 

B ( J ) =- (DYBU1 ♦ VV* ( BX 1 <I)>BX2(I) > ) 

C ( J) =0. 

D(J)=DYBU2*FPU(I)-VV*(BX1 (I)*PHI <LR> ♦BX2(I)*PHI (LL) > 

GO  TO  200 
260  CONTINUE 

IF  < V ( J) .GT.O.)  GO  TO  265 
C HYPERBOLIC 

OMEGA (J)=OMEGAH 

VV=K-CX(I-1)*(PHI (L)-PHI (LL) ) -CX ( 1-2) * (PHI (LL)-PHI (LLL) ) 

A ( J ) =DYBU1 

B ( J ) =VV*BX 1 (I-l)-DYBUl 
C ( J) =0. 

D ( J)=DYBU2*FPU( I ) ♦VV* (BX1 ( 1-1 ) *PHI (LL) *BX2( 1-1 ) * 

1 (PHI (LL)-PHI (LLL) ) > 

GO  TO  200 
C PARABOLIC 
265  CONTINUE 

OMEGA ( J)=OMEGAP 
A ( J ) =DYBU1 

B ( J ) =VV*BX 1 ( 1-1 )-DYBUl 
C(J)=0. 

D ( J) =DYBU2*FPU< I ) ♦VV* (BX1 ( 1-1 )*PHI (LL) *BX2 ( 1-1 ) * 

1 (PHI (LL)-PHI (LLL) ) ) 

GO  TO  200 

C BODY  BOUNDARY  BELOW 
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270  CONTINUE 

VV=K-AX1 (I)*(PHI (LR)-PHI (L) ) -AX2 ( I ) * (PHI (L)-PHI <LL> ) 

IF  (VV.LT.O.)  GO  TO  280 
C ELLIPTIC 

OMEGA (J)=OMEGAE 
A(J)=0. 

B( J)=-(DYBL1*VV*(8X1 <I)*BX2(I>  ) ) 

C( J)=0YBL1 

D(J)=-0YBL2*FPL(I)-VV*(BX1 (I)*PHI <LR> ♦BX2(I)*PHI (LL) ) 

GO  TO  200 
280  CONTINUE 

IF  (V(J) .GT.O.)  GO  TO  285 
C HYPEPBOLIC 

OMEGA (J)=OMEGAH 

VV=K-CX(I-I)*(PHI (L) -PHI (LL ) ) -CX ( 1-2 ) * (PHI (LL)-PHI (LLL) ) 
A(J)=0. 

B(J)=VV*BXl (I-l)-DYBLl 
C( J)=DYBL1 

D ( J>  =-DYBL2*FPL ( I ) ♦VV* (BX1 (I-I)*PHI (LL > *BX2 ( I- 1 > * (PHI (LL)- 
1 PHI (LLL))) 

GO  TO  200 
C PARABOLIC 
285  CONTINUE 

OMEGA ( J) =OMEGAP 
A ( J) =0 • 

B(J)=VV*BXI (I-l)-DYBLl 
C( J)=0YBL1 

D( J)=-DYBL2*FPL(I > ♦VV* (BX1 ( 1-1 >*PHI ( LL ) ^BX2 ( I - 1 ) * 

1 (PHI (LL)-PHI (LLL) ) ) 

GO  TO  200 
290  CONTINUE 
C BODY  BOUNDARY  J=JW 
A ( J) =0* 

B ( J) =1  . 

C ( J) =0. 

D ( J) =PHI (L) 

200  CONTINUE 

PHI (LR) =TPHIR 
PHI (LL ) =TPHIL 
160  CONTINUE 

C TRIDIAGONAL  MATRIX  IS  SET  NOW  SOLVE  FOR  COLUMN  OF  PHI 
CALL  TRI  (I) 

C CHECK  FOR  COLUMN  CONVERGENCE  OF  PHI 
DO  295  J=2* JM 1 
L=M*  J 
JERROR=J 

ERRC=PHIOC ( J) -PHI (L) 

IF  (ABS(ERRC) .GT.EPSCOL)  GO  TO  150 
295  CONTINUE 
294  CONTINUE 

IF  (NITERC.GT.NCOL)  WRITE  (6»904)  NI TERG* I *ERRC» JERROR 
C CONVERGED^  RELAX  PHI,  FIND  ERROR,  CAL.  V,  AND  MOVE  TO  NEXT  COLUMN 
DO  296  J=2 ♦ JM 1 
L=M  ♦ J 

ERR = OMEGA ( J) * (PHI (L)-PHIOG(J) ) 
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PHI ( L ) =PH I OG ( J ) ♦ERR 

V(J)=K-AX1 ( I ) * ( PHI (L^JM)-PHI <L > > -AX2 ( I > * (PHI (L > -PH I (L-JM) ) 
IF  (ABS(ERR) .LT.ARS(ERROR) ) GO  TO  296 
ERROR=ERR 
LERROR=L 
296  CONTINUE 

IF  (IFLAG.NE.l)  Go  TO  100 
L=M*JW 

PHI (L)=PHI (L-l ) ♦OY( JWM1)*(PHI (L-l)-PHI (L-2) ) /OY ( JW-2) 

PH I UB ( I ) =PH I (L>l)-OY(JW)*(PHI (L^2)-PHI (L> 1 ) ) /DY ( JWP1 ) 
V(JW)=K-AX1 (I)* (PHI (L,JM)-PHI (L) ) -AX2 ( I ) * (PHI (L)-PHI (L-JM) ) 
IF  (I.EQ.ITE)  GAMTE=PHIUB ( I ) -PHI (L ) 

IF  (ABS(GAMFF) .LE.l.E-8)  GAMTE=0. 

IF  (K.LT.O.)  GAMFF=GAHTE 

IF  (I.EQ.ITE. AND. NITERG.E0.1)  GAMTE1=GAMTE 
100  CONTINUE 

PRINT  OUT  ERROR  AFTER  GRID  SWEEP 

WRITE  (6*905)  NITERG.ERROR.LERROR 
IDOU8=0 

IF  ( ABS (ERROR) .LE.EPSGRD (KGRD) ) GO  TO  300 
IF  (NITERG.EQ.NGRID)  GO  TO  310 
IF  (MOD (NITERG, NDUMP) .EQ.O)  GO  TO  310 
GO  TO  50 

300  CONTINUE 
KGRD=KGRD* 1 
IDOUB= 1 

GO  TO  310 

301  CONTINUE 

IF  (K.LT.O.) 

CALL  GAMFUN 
WRITE  (6.910) 

302  CONTINUE 
CALL  FPRINT 
WRITE  (6,900) 

WRITE  (6,906) 

CALL  DOUBLE 
WRITE  (6,914) 

(6.902) 

(6.903) 

(6,911) 

(6,903) 

50 


GO  TO  302 

niterg,gamff,gamte,doublt 


niterg 


(X(l) ,I*1,IM) 


(Y(I) ,I=1,JM) 


IM,JM,JW,ILE»ITE 

WRITE 
WRITE 
WRITE 
WRITE 
GO  TO 
310  CONTINUE 
C TAPE  DUMP 

WRITE  (7)  NITERG, I M ♦ I M 1 ♦ JM » JM1 , JW , JWP 1 » JWM 1 , I TE » ILE »K , DEL, 

1 ALPHA, GAMTE1 »GAMFF ,NDB , XH ♦ M8.GAM ,DYBU1 ,DYBU2,DYBL1,DYBL2, 

2 DOUB.DOUBLT.ALPHAF 


WRITE 
WRITE 
WRITE 
WRITE 
WRITE 
L=IM*JM 
WRITE  (7) 
END  FILE 


(7) 

(7) 

(7) 

(7) 

(7) 


(X ( I ) ♦ 1=1 ,IM) 

< Y ( I ) .1  = 1 . JM) 

( FPU ( I ) ,FPL(I) »PHIUB(I) *I=ILE*ITE) 

(AX1 (I) ,AX2(I) ,BX1 (I) , BX2 ( I ) ,CX(I) ,DX(I)  ,1  = 1 ,IM1 ) 
( AY  1 ( I ) , AY2 ( I ) ,DY(I) ,I  = 1,JM1) 

(PHI (I) ,1=1, L) 
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WRITE  (6*907)  NITERG 
CALL  PRINT (NITERG) 

IF  (KGRD.GT.KEPS)  GO  TO  320 
IF  (NITERG. EQ.NGRID)  GO  TO  330 
IF  ( IDOUB.EQ. 1 ) GO  TO  301 
GO  TO  50 
320  CONTINUE 

WRITE  (6*908)  NITERG 
GO  TO  350 
330  CONTINUE 

WRITE  (6*909)  NITERG 

900  FORMAT  (1H1) 

901  FORMAT  (1H  ♦/**  CASE  IS  BEING  RESTARTED  AT  I TERAT ION* 15 ) 

902  FORMAT  (1H  */,*  X<I) ,I=1*IM*) 

903  FORMAT  (10E13.5) 

904  FORMAT  (1H  */**  AT  I TERAT I ON* 15*  C0LUMN*I4*  FAILED  TO  CONVERGE* 

1 * ERR  =*E 13.5*  J=* 1 3 ) 

905  FORMAT  (1H  */♦*  AT  I TERAT ION* 1 5*  THE  MAXIMUM  ERROR  =*E13.5 

1 * AND  OCCURRED  AT  N0DE*I5) 

906  format  ( i h */**  the  number  of  nodes  is  being  doubled  at  iteration* 

1 15) 

907  FORMAT  (1H  */,*  TAPE  HAS  BEEN  DUMPED  AT  I TERAT I ON* 1 5 > 

908  FORMAT  (1H  ,/,*  SOLUTION  HAS  CONVERGED  TO  DESIRED  ACCURACY  AT  ITER 
1 AT I0N*I5) 

909  FORMAT  (1H  ,/**  MAXIMUM  NUMBER  OF  ITERATIONS  HAS  BEEN  REACHED,  CAS 
IE  IS  BEING  TERMINATED  AT  ITERATI0N*I5) 

910  FORMAT  (1H  ,/,*  UPDATE  GAMFF  AND  FARE IELD  AT  ITERATION*I5 

1 * GAMFF  =*E 13.5*  GAMTE  =*E13.5*  DOUBLET  =*E13.5) 

911  FORMAT  ( 1H  */,*  Y ( I ) * 1 = 1 * JM*) 

912  FORMAT  (8A10) 

913  FORMAT  (1H  ,/**  SIMILARITY  PARAMETER  (K)  =*E13.5,/,*  SCALING  FACTO 
1R  ( CP/CPBAR ) =*E 13.5) 

914  FORMAT  (1H  */**  IM  =*I4*  JM  = *I4*  JW  = *I4*  ILE  =»I4*  ITE  =*I4) 

350  CONTINUE 

CALL  FPRINT 
END 

SUBROUTINE  DOUBLE 
REAL  K * M8 

COMMON  /DELTA/  DX ( 99 ) * DY ( 99) * AX  1 (99) * AX2 ( 99 ) *BX 1 ( 99) * BX2 (99) * 

1 CX (99) , AY1 (99) ♦ AY2 (99) *X ( 100) *Y ( 100) * FPU (99) * FPL (99) ,PHIUB (99) * 

2 IM* IM1 * JM, JM1 * JW* JWP1 * JWM1 ,ITE, ILE*DYBU1 *DY0U2»DYBL1 *DYBL2» 

3 ALPHA *DEL»M8*GAM*K*NDB*XH» TITLE (8) » D0U8  * DOUBLT  * ALPHAF 
COMMON  /COEFF/  A (99) *B (99) ,C (99) *D (99) *PHI ( 10000) 

COMMON  /GAMMA/  GAmTEI , GAMTE *PGFF *GAMFF 

NDB=NDB* 1 

C CHANGE  GRID  INDICES 
I MAX=2* I M-4 
I MAX  1 = I MAX- 1 
JMAX=2* JM-3 
JM AX  1 = JMAX  — 1 
JWN=2* JW-2 
JWNP 1 = JWN* 1 
JWNM 1 = JWN- 1 
ITEN=2*I TE-3 

D I ST= I •/ ( 1 .♦DY (JWPl ) /DY (JW) ) 
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IK= I TE 
C CHANGE  X 
L = IM 

DO  30  1=5* IMAXl *2 
II=IMAX1-I*5 
L=L- 1 

X ( 1 1 ) =X (L ) 

X(II-1)*.5*(X(L)«X(L-1)) 

30  CONTINUE 

X(IMAX)=2.*X<IMAXl)-X UMAX  1-1) 

DELT=X(4)-X(3> 

X ( 2 ) =X ( 3) -DELT 
X ( 1 ) =X (2) -DELT 
C FIND  NEW  LEADING  EDGE 
DO  31  1=3* IMAX 1 
IF  ( X ( I ) .GE.O.)  GO  TO  32 

31  CONTINUE 

32  ILENM 
C CHANGE  Y 

L = JM 

DO  40  1=2* JHAXl *2 
1 I = JMAX 1 -1*2 
L=L- 1 

Y ( 1 1 ) =Y (L> 

Y(II-1)=.5*(Y(L)*Y(L-1>  ) 

IF  (II-1.EQ.JWNP1)  Y ( 1 I - 1 ) =Y  < L — 1 ) ♦DIST* ( Y (L ) - Y (L-l ) ) 

IF  (II-1.EQ.JWNM1)  Y<II-1)=Y(L)*DIST*(Y<L-1)-Y<L) ) 

40  CONTINUE 

Y ( JMAX >=2.*Y( JMAX 1 >-Y( JMAX 1-1 ) 

Y< 1)=2.*Y(2)-Y(3) 

C MOVE  THE  PHI (S)  (INTERIOR  ONLY)  AND  INTERPOLATE  FOR  PHI  <I=ODD> 
L=IM1* JM 

DO  10  1=3* IMAXl *2 
II=IMAX1-I*3 
LL= ( 1 1-1 ) * JMAX 
DO  15  J=2  * JMAX 1*2 
JJ= JMAX 1-J»2 
LLL=LL4 JJ 
L=L- 1 

PHI (LLL ) =PHI (L) 

PHI (LLL-1)=.5* (PHI (L) ♦PHI (L-I ) ) 

IF  (JJ.NE.JWN)  GO  TO  15 

PHI (LLL-l ) =PHI (L> ♦DIST* (PHI (L-l) -PHI  (L) ) 

IF  (II.GT.ITEN)  GO  TO  16 
IF  (II.GE.ILEN.AND.II.LE.ITEN)  GO  TO  17 
PHI (LLL* 1 ) =PHI (L ) *DIST* (PHI (L*1)-PHI (L) ) 

GO  TO  15 

16  CONTINUE 

SIGI=(X(II)-1.)*(GAMFF-GAMTE)/(X(IMAX1)-1.)*GAMTE 
PDUM=PHI (L) *.5*5101 

PHI (LLL*1 ) =PDUM*DIST* (PHI (L*1)-PDUM) 

PDUM=PH I (L)-.5*SIGI 

PHI (LLL-l ) =PDUM*DIST* (PHI (L-l)-PDUM) 

GO  TO  15 

17  CONTINUE 
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PHI  (LLL*1)=PHIUB(IK)  ♦DISTMPHI  < L ♦ 1 > -PHIUB  ( IK ) ) 

I  K= I K- 1 
15  CONTINUE 
L=L-2 

10  CONTINUE 

C FILL  IN  OTHER  PHI(S)  (I=EVEN) 

IMAX2= I MAX  1-1 
DO  20  I=4» IMAX2*2 
I I = IMAX2-I*4 
L= ( 1 1 - 1 ) * JMAX 
DO  25  J=2 ♦ JMAX 1 
JJ= JMAX 1 -J*2 
LL=L*JJ 

PHI (LL)=.5*(PHI (LLOMAX)  *PHI (LL- JMAX ) ) 

25  CONTINUE 
20  CONTINUE 
IM=IMAX 
IM 1 = I MAX  1 
JM=JMAX 
JM 1 = JM  AX  1 
JW= JWN 
JWP 1 = JWNP1 
JWM 1 = JWNM 1 
ITE=ITEN 
ILE= ILEN 

KTE=(ITE-1)*JM*JW 

PHI (KTE*JM)=.5*PHI (KTE) ♦.25*GAMTE*.5*PHI (KTE*2*JM> 

C INITIALIZE  FINITE  DIFFERENCE  COEFFICIENTS  AND  F ARFIELD 
CALL  INITAL 
CALL  FARFLD 
L=(ILE-1)*JM*JW 

PHIUB ( ILEN)=PHI (L ♦ 1 ) ~DY ( JW ) * ( PHI <L*2)-PHI ( L ♦ 1 > ) /DY ( JWP1 ) 

RETURN 

END 

SUBROUTINE  FARFLD 
REAL  K »M0 

COMMON  /DELTA/  DX ( 99 ) *DY (99) » AX  1 (99) » AX2 ( 99) ♦ BX1 (99) *BX2 ( 99 ) ♦ 

1 CX(99) »AY1 (99) *AY2(99) * X ( 1 00 ) * Y ( 1 00 ) *FPU(99> .FPL (99) .PHIUB (99) » 

2 IM» IM1 » JM* JM1 .JW.JWP1.JWM1.ITE* ILE.DYBU1 .DYBU2.DYBL1 .DYBL2* 

3 ALPHA.DEL.M8.GAM.K.NDB.XH.TITLE (8) .DOUB.DOUBLT . ALPHAF 
COMMON  /COEFF/  A(99) *B (99) .C (99) *D (99) .PHI (10000) 

COMMON  /GAMMA/  GAMTE 1 ♦ GAMTE . PGFF ♦ GAMFF 

IF  (K.LT.O.)  GO  TO  25 
C SUBSONIC  FARFIELD 
SUM=0 • 

DO  1 I ♦ I M 1 
SUM  1 =0  • 

M= ( I - 1 ) * JM 

PO=PHI (M*2)-PHI <M-JM*2) 

IFLAG=0 

IF  ( I .GE  » I LE» 1 .AND. I.LE.ITE)  IFLAG=1 

DO  2 J=3  » JM1 
L = M ♦ J 

IF  (IFLAG.EQ.1.AND.J.EQ.JWP1)  PO=PHlUB ( I ) -PHIUB ( 1-1 ) 

PN=PHI (L)-PHI (L-JM) 
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IF  (J.EQ.JW.AND.I.EQ.ITEM)  PN=PN- »5*GAMTE 1 
SUM1=SUM1 ♦ (PN*PO) **2*DY (J-l ) 

PO=PN 

2 CONTINUE 

SUM = SUM ♦ «25*SUMl/DX  ( 1-1 ) 

1 CONTINUE 
SUM=.25*SUM 
DOUBLT=DOUB*SUM 
SK=SQRT (K ) 

C0N=D0UBLT/(3.14159265*SK) 

CONl=l ./6. 2831853 
C0N2=1. 570796325 
DO  10  1=1 » IM 
M=(I-1)*JM 
J3= JM 1 

IF  (I.LE.2.0R.I.E0.IM)  J3=l 
DO  15  J=1 » JMt J3 
L=M>  J 

IF  ( Y ( J ) .EQ.O.)  GO  TO  20 

PHI (L ) =CON*X ( I ) / (X ( I > **2*K*Y ( J) **2) ♦GAMFF*C0N1« ( AT AN (X ( I > 

1 / ( SK*Y (J) ) ) ♦SIGN (C0N2  * Y ( J) ) ) 

GO  TO  15 
20  CONTINUE 

PHI ( L ) =CON/X  ( I ) 

15  CONTINUE 
10  CONTINUE 
RETURN 
25  CONTINUE 

C SUPERSONIC  FARFIELD 
DO  30  1 = 1 * IM 
M= ( I - 1 ) * JM 
J3= JM 1 

IF  (I.LE.2.0R.I.E0.IM)  J3=l 
DO  35  J= 1 * JM  * J3 
L=M  ♦ J 
PHI ( L ) =0 • 

35  CONTINUE 
30  CONTINUE 
RETURN 
END 

FUNCTION  FL  (CAC) 

REAL  K *M8 

COMMON  /DELTA/  DX(99) *DY (99) »AX1 (99) *AX2(99) *BXl (99) *8X2(99) * 

1 CX (99) ♦ AY1 (99) *AY2 (99) *X ( 100) *Y < 100) *FPU (99) *FPL(99) *PHIUB(99) t 

2 IM»IM1»JM*JM1 * JWt JWP1 ♦ JWM1 * ITE* ILEtDYBUl ♦DYBU2*DYBL1 *DYBL2» 

3 ALPHA »DEL»M8*GAM*K»NDB*XH* TITLE (8) *DOUB*DOUBLT » ALPHAF 
C AIRFOIL  LOWER  SURFACE  SHAPE  FUNCTION 

FL=-.5>2.* ( CAC-. 5 )**2-ALPHA«C AC/DEL 

IF  (CAC.GT.XH)  FL=FL~ALPHAF  * (CAC-XH) /DEL 

RETURN 

END 

FUNCTION  FLP  (CAC) 

REAL  K »M8 

COMMON  /DELTA/  DX ( 99) *DY ( 99 ) ♦ AX1 ( 99) ♦ AX2 ( 99) ♦ BX 1 ( 99 ) »BX2 (99 ) ♦ 

1 CX(99) * A Y 1 (99) ♦ AY2 (99 ) * X ( 1 00 ) *Y (100) * FPU (99) * FPL (99) *PHIU8(99) t 
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2 IMtlMl t JM* JM1 » JW* JWP1 * JWM1 * ITE* ILE *DYBU1 *DYBU2 *DYBL  1 *DYBL2» 

3 ALPHA *DEL.M8*GAM*K*NDB*XH* TITLE (8) *D0UB*D0UBLT  *ALPHAF 
C AIRFOIL  LOWER  SURFACE  SLOPE  DISTRIBUTION 

FLP=A.*(CAC-. 5) -ALPHA/DEL 

IF  (CAC.GT.XH)  FLP=FLP-ALPHAF /DEL 

RETURN 

END 

SUBROUTINE  FPRINT 
REAL  K »M8 

COMMON  /DELTA/  DX ( 99) ♦ DY ( 99 > ♦ AX  1 ( 99 > * AX2 < 99) ♦ BX 1 ( 99 > * BX2 <99 > , 

1 CX<99) »AY1 (99) *AY2<99) *X<100) * Y < 1 00 ) *FPU < 99) *FPL<99) *PHIUB<99) ♦ 

2 I M ♦ I M 1 * JM  ♦ JM  1 *JW*JWP1 » JWM 1 * ITE* ILE*DYBU1 ♦ DYBU2 * DYBL 1 * DYBL2  * 

3 ALPHA ♦DEL*M8»GAM*KtND8»XHt TITLE (8) *DOUB*DOUBLT  * ALPHAF 
COMMON  /COE FF/  A (99) *B (99) *C (99) *D (99) ,PHI ( 10000) 

COMMON  /GAMMA/  GAMTE 1 * GAMTE * PGFF * GAMFF 

CPCPB=DEL**. 6666666667/ < ( 1 . ♦GAM ) *M8**2> **. 3333333333 

PART=.5*< (X(ITE*1)-1.)*(GAMFF-GAMTE)/(X(IM1)-1.)4GAMTE) 

I JW=ITE*JM*JW 

PHI (IJW)=PHI (I JW ) -PART 
L=(ILE-2)*JM*JW 
PHIUB(ILE-1)=PHI (L) 

PHIUB  ( ITEM  >=PHI  ( I JW)  *2. "PART 
C COMPUTE  CP  LOWER  (A)  AND  CP  UPPER  (B) 

DO  10  I = I LE  * I TE 
M= ( I — 1 ) * JM 
L=M* JW 

A ( I ) = -2 • * ( AX  1 ( I ) * ( PH  I (L»JM)-PHI (L ) ) ♦ AX2 ( I ) * (PH  I (L)-PHI (L-JM) ) ) 

1 *CPCPB 

B ( I ) =-2.* (AX1 ( I ) * (PHIUB ( I ♦ 1 ) -PHIUB ( I ) ) ♦AX2(I)*( PHIUB (I)- 
1 PHIUB(I-I) ) )*CPCPB 
10  CONTINUE 

PHI <IJW)=PHI ( I JW) *PART 
C COMPUTE  AIRFOIL  FORCE  COEFFICIENTS 
C10=A(ILE)-B(ILE) 

CLIFT=C10*X(ILE) 

C20=CL I FT 

CM0M=.5*CLIFT*X(I(  E) 

C3O=0 • 

CHINGE=C30 

IF  (XH.GE.X(ILE) ) GO  TO  25 
C30=C10*(X(ILE)-XH) 

CHI NGE= • 5*C30* ( X ( ILE ) “XH ) 

25  CONTINUE 
ILE1=ILE*1 
DO  30  I=ILE1 ♦ ITE 
C 1 =A ( I ) -B ( I ) 

C2=C 1 *X ( I ) 

IF  (X(I).GT.XH)  C3=C1*(X(I)-XH) 

CLIFT =CL IFT ♦•5*(C1*C10)*DX(I-1) 

CMOM=CMOM> .5* (C2*C20) *DX  ( 1-1 ) 

DXX=DX ( 1-1 ) 

IF  ( X ( I ) .GT.XH.AND.X(I-l) .LE.XH)  DXX=X(I)-XH 

CHINGE=CHINGE*.5* (C3*C30)*DXX 

C10=C1 

C20=C2 


-103- 


30 


C30=C3 

CONTINUE 

CIRLIF=2.*GAMTE*CPCPB 

CPCRIT=-2.*CPCPB*K 


WRITE 

(6*900) 

WRITE 

(6*901) 

(TITLE(I) *1 

00 

it 

WRITE 

(6*902) 

M8 

WRITE 

(6*903) 

K 

WRITE 

(6,904) 

DEL 

WRITE 

(6*905) 

ALPHA 

WRITE 

(6*916) 

XH 

WRITE 

(6,917) 

ALPHAF 

WRITE 

(6,906) 

CPCPB 

WRITE 

(6,918) 

CPCRIT 

write 

(6*907) 

gamte 

WRITE 

(6*908) 

cirlif 

write 

(6*909) 

WRITE 

(6*910) 

CLIFT  *CMOM  * 

CHINGE 

WRITE 

(6*911) 

WRITE 

(6,912) 

WRITE 

(6,915) 

(X  < I ) » 1 = ILE 

* I TE ) 

WRITE 

(6*913) 

WRITE 

(6,915) 

(B  ( I ) *1  = ILE 

* I TE ) 

WRITE 

(6*914) 

WRITE 

(6*915) 

( A ( I ) * I = ILE 

, I TE ) 

FORMAT 

( 1 H 1 ) 

FORMAT 

( 30X  * 8A 1 0 ) 

FORMAT 

< 1 H ,/ 

♦1H  * / ♦ 1 H »/ 

,*  MAC 

FORMAT 

(*  SIM 

ILARITY  PARAMETER 

FORMAT 

(*  THI 

CKNESS  RATIO 

= *E  1 3 

900 

901 

902 

903 

904 

905  FORMAT  (* 

906  FORMAT  (* 

907  FORMAT  (* 

908  FORMAT  (* 
1RCULATION) 

909  FORMAT  (1H 

910  FORMAT  (1H 
1 3X*H I NGE 

911  FORMAT  <1H 
1ED ) * ) 

912  FORMAT 

913  FORMAT 

914  FORMAT 

915  FORMAT 

916  FORMAT 

917  FORMAT 

918  FORMAT 
RETURN 
ENO 

FUNCTION  FU 


H NUMBER  =*E13«5) 
(K)  =*E13.5> 


AIRFOIL  ANGLE  OF  ATTACK  (RAOIANS)  =*E13.5> 

CP  SCALING  FACTOR  (CP/CPBAR)  =*E13.5> 

SCALED  AIRFOIL  CIRCULATION  =*E13.5) 

LIFT  COEFFICIENT  BASED  ON  C IRCULAT ION*27H 
=E13.5) 

* / * 1 H ,/**  AIRFOIL  FORCE  COEFFICIENTS*) 

♦ / ♦ 3X*L IFT  =*E13.5,/*3X*M0MENT  ABOUT  (X=0)  =*E13.5*/* 
MOMENT  =*E13.5) 

• / * 1 H ,/,*  PRESSURE  COEFFICIENTS  ON  THE  AIRFOIL  (UNSCAL 


<2*CP/CPBAR*CI 


(1H 

(1H 

(1H 


*/*3X*AlRFOIL  COORDINATE  =*) 

» / * 3X*A IRFO IL  PRESSURE  COEFFICIENTS* 
*/*3X*AlRF0IL  PRESSURE  COEFFICIENTS* 
(3X10E13.5) 

(*  HINGE  POINT  =*EI3.5) 

(*  FLAP  ANGLE  (RADIANS)  =*E13.5) 

(*  CRITICAL  PRESSURE  COEFFICIENT  (SONIC) 


(CAC) 


UPPER 

LOWER 


= *) 
=*) 


=*E13.5) 


REAL  K ♦ M8 

COMMON  /DELTA/  DX (99) *DY (99) ♦ AX1 (99) ♦ AX2 (99) *BX1 (99) *BX2 (99) , 

1 CX (99) * AY1 (99) *AY2(99) * X ( 1 00 ) ♦ Y ( 1 00 ) *FPU ( 99) *FPL(99) *PHIUB(99) ♦ 

2 IM*IM1 ♦ JM*JM1 * JW*JWP1 ♦ JWM1 , I TE * ILE *DYBU1 ♦ DYBU2 ♦ DYBL 1 *DYBL2* 

3 ALPHA ♦ DEL *M8*GAM*K»NDB*XH* TITLE (8) *DOUB*DOUBLT  *ALPHAF 
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C AIRFOIL  UPPER  SURFACE  SHAPE  FUNCTION 
FU=.5-2.*(CAC-.5)**2-ALPHA*CAC/DEL 
IF  (CAC.GT.XH)  FU=FU-ALPHAF*(CAC-XH)/OEL 
RETURN 
END 

FUNCTION  FUP  (CAC) 

REAL  K*M8 

COMMON  /DELTA/  DX ( 99 > *DY (99) * AX  1 (99) * AX2 (99) • BX 1 (99 > ♦ BX2 ( 99 ) . 

1 CX(99) ♦ AY  1 (99) »AY2(99) ♦ X ( 1 00 ) ♦ Y (1 00 ) *FPU < 99) ,FPL (99) *PHIUB(99) » 

2 IM. IM1 * JM* JM1 . Jw* JWP1 » JWM1 ♦ ITE* ILE*DYBU1 *DYBU2*DYBL1 .DYBL2* 

3 ALPHA  t DEL *M8*GAM*K  *NDB *XH*TITLE(8) *DOUB*DOUBLT * ALPHAF 
C AIRFOIL  UPPER  SURFACE  SLOPE  DISTRIBUTION 

FUP=-4.*(CAC-.5> -ALPHA/DEL 
IF  (CAC.GT.XH)  FUP=FUP-ALPHAF /DEL 
C DOUB  IS  DOUBLET  STRENGTH  DUE  TO  THICKNESS 
DOUB=. 6666666667-2. *DEL**2 
DOUBLT=DOUB 
RETURN 
END 

SUBROUTINE  GAMFUN 

COMMON  /GAMMA/  GAMTE 1 tGAMTE tPGFF * GAMFF 
GAMFF=GAMTE1*PGFF* (GAMTE-GAMTE1 ) 

GAMTE 1 =GAMTE 

RETURN 

END 

SUBROUTINE  INITAL 
REAL  K.M8 

COMMON  /DELTA/  DX ( 99 > . DY ( 99 > . AX  1 ( 99) ♦ AX2 < 99 ) * BX 1 ( 99 ) ♦ BX2 ( 99 > , 

1 CX(99) »AY1 (99) ,AY2(99) tX(lOO) »Y(100) » FPU (99) , FPL (99) .PHIUB(99) . 

2 IM. IM1 * JM* JM1 ♦ JW* JWP1 ♦ JWM1 ♦ ITE* ILE*DYBU1 *DYBU2*DYBL1 .DYBL2. 

3 ALPHA ♦ DEL ♦ MR ♦ GAM *K«NDB»XH* TITLE (8) * DOUB tDOUBLT* ALPHAF 
C CALCULATE  DX 

DO  15  I = 1 t I M 1 
DX ( I ) =X ( I ♦ 1 ) -X ( I ) 

15  CONTINUE 
C CALCULATE  DY 

DO  25  1 = 1*  JM 1 
DYII)*Y(I*1)-Y(I) 

25  CONTINUE 

DO  30  I =2 ♦ I M 1 

AX  1 (I)=DX(I-1)/(DX(I)*(DX(I-1)*DX(I) ) ) 
AX2(I)=DX(I)/(DX(I-1)*(DX(I-1) ♦DX(I) ) ) 

BX 1 ( I ) =2.*AX1 ( I ) /DX ( I - 1 ) 

BX2(I)=2.*AX2(I)/DX(I) 

CX ( I ) = ,5/DX ( I ) 

30  CONTINUE 

CX ( 1 ) =.5/DX ( 1 ) 

DO  40  I =2 1 JM1 

AY1 (I)=2./(DY(I)*(DY(I)*DY(I-1) )) 

AY2(I)=2./(DY(I-1)*(DY(I)^DY(I-1) ) ) 

40  CONTINUE 

IF  (K.GT.O.)  GO  TO  50 

C TOP  AND  BOTTOM  AND  RHS  BOUNDARY  CONDITIONS 
AX  1 ( IM1 ) =0. 

AX2 ( IM1 ) =0. 
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BX1 ( IM 1 ) =0 • 

BX2(IM1)=2./DX(IM1-1>**2 
AYl <2)=2./DY<2>**2 
A Y2 ( 2 ) =0 • 

AY  1 ( JM 1 ) =0 • 

AY2 ( JMl ) =2./DY ( JM1-1 ) **2 
50  CONTINUE 

DYBUl=2./< (DY(JWPl) ♦2.*0Y(JW) )*DY( JWP1) ) 

DYBU2=DY( JWP1 >*DYBUl 

DYBLl=2./<  <DY(JW-2) ♦2.*DY(JWM1 ) >*DY( JW-2) ) 

DYBL2=DY ( JW-2) *DYBL1 
C SET  AIRFOIL  BOUNDARY  CONDITION 
DO  45  I=ILE* ITE 
FPU ( I ) =FUP ( X ( I ) ) 

FPL<I)=FLP(X(I) ) 

45  CONTINUE 
RETURN 
END 

SUBROUTINE  PRINT  (NITERG) 

REAL  K,M8 

COMMON  /DELTA/  DX (99) ,DY (99) ♦ AX1 (99) ,AX2 (99) »BXl (99) ,BX2 (99) * 

1 CX(99) »AY1 (99) ,AY2(99) * X ( 1 00 ) ♦ Y ( 1 00 ) »FPU(99) , FPL (99) ,PHIUB(99) , 

2 IM, IM1,JM, JM1 »JW, JWP1, JWMl ,ITE,ILE»DYBU1»DYBU2,DY8L1 ,DYBL2, 

3 ALPHA »DELtM8.GAM.K»NDB»XH» TITLE (8) »DOUB»DOUBLT ♦ ALPHAF 
COMMON  /COEFF/  A(99) *B ( 99 ) »C (99 ) *D ( 99) »PHI (10000) 

COMMON  /GAMMA/  GAmTEI »GAMTE  tPGFF  »GAMFF 

PART  = .5*(  (X(ITEM)-1.)*(GAMFF-GAMTE)/(X(IM1)-1.)*GAMTE> 

I JW=ITE*JM> JW 
PHI ( I JW) =PHI ( I JW) -PART 
L= ( ILE-2) * JM* JW 
PHIUB(ILE-l) =PHI (L) 

PHIUB(ITE*1)=PHI ( I JW) *2.*PART 
C COMPUTE  CP  LOWER  (A)  AND  CP  UPPER  (B) 

DO  297  I*  ILE  * I TE 
M= ( I - 1 ) * JM 
L=M ♦ JW 

A ( I ) = -2 • * ( AX  1 (I)* (PHI (L* JM) -PHI (L ) ) ♦ AX2 ( I ) * ( PH  I (L)-PHI (L-JM) ) ) 

B ( I ) = -2 • * ( AX  1 (I)*(PHIUB(I*1)-PHIUB<I> ) *AX2 ( I ) * (PHIUB ( I ) - 
1 PHIUB ( 1-1 ) ) ) 

297  CONTINUE 

PHI ( I JW) =PHI ( I JW) *PART 

WRITE  ( 6 » 9 1 1 ) NITERG 

WRITE  (6,903)  ( B ( I ) , 1= ILE ♦ I TE ) 

WRITE  (6,912)  NITERG 

WRITE  (6,903)  ( A ( I > , I = ILE , I TE > 

903  FORMAT  (10E13.5) 

911  FORMAT  (1H  ,/,*  AT  I TERAT ION* 15*  SCALED  PRESSURE  COEFFICIENT,  UPPE 
1R  (ILE  TO  ITE)  =*) 

912  FORMAT  (1H  ,/,*  AT  I TERAT I ON* 15*  SCALED  PRESSURE  COEFFICIENT,  LOWE 
1R  (ILE  TO  ITE)  =*> 

RETURN 

END 

SUBROUTINE  TRI  (I) 

REAL  K»M8 

COMMON  /DELTA/  DX (99) ,DY (99) , AXl (99) ,AX2 (99) ,BX1 (99) ,BX2 (99) , 
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1 CX  ( 99 ) ♦ AY  1 (99) *AY2  C 99 > ♦ X ( 1 00 > » Y < 1 00 ) *FPU<99>  tFPL(99) ,PHIUB(99) t 

2 IM» IM1 ♦ JM* JM1 ♦ JWt JWP1 ♦ JWM1 » ITE*ILE*DYBU1 ♦DYBU2.DYBL1 .DYBL2* 

3 ALPHA ♦0ELtH8*GAM»K*NDB.XHt TITLE (8) *DOUB»DOUBLT  * ALPHAF 
COMMON  /COE FF/  A (99) »B (99) ,C (99) .0 (99) ,PHI ( 10000) 

DO  10  KK=3 ♦ JM 1 

J=JM1-KK*3 

P=A(J-1)/B(J> 

B(J-1)=B(J-1)-P*C(J> 

D ( J-l ) =D ( J-l ) -P*D ( J) 

10  CONTINUE 
M= ( I- 1 ) * JM 
PHI (M*2)=D(2)/B(2) 

DO  20  J=3? JM1 
L=M*J 

PHI ( L ) = ( D ( J) -PHI (L-1)*C( J) )/B(J) 

20  CONTINUE 
RETURN 
END 
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APPENDIX  B 


FORTRAN  LISTING  OF  UTRANS 


A FORTRAN  listing  of  the  source  deck  for  the  UTRANS 
program  is  presented  in  the  following  pages.  The  program, 
as  configured  here,  requires  144 8K  words  to  load  and  132 0K 
words  to  execute. 


-108- 


PROGRAM  UTRANS  < INPUT*0UTPUT*TAPE5=INPUT*TAPE6=0UTPUT,TAPE7* 

1  TAPE8*TAPE9) 

COMPLEX  PHIUB*B*D,PHI,GAMTEl*GAMTE*GAMFF.PHIOG*PART.ERR, 

1 ERR0R*0MEG2I*TPHIL*TPHIR,SIGI*GAMFFS*C0NAFF,PHIAFF 
REAL  K *M8 

COMMON  /DELTA/  DX ( 99) *DY ( 99) ♦ AX  1 (99) ♦ AX2 < 99) *BX 1 (99) *BX2 (99) « 

1 CX(99) *AY1 (99) *AY2(99) »X<100) • Y ( 1 00 ) • FPU ( 99) *FPL<99) *PHIUB«99) * 

2 IM»IM1*JM* JM1 * JW* JWP1* JWM1*ILE*ITE*DYBU1*DYBU2*DYBL1*0YBL2* 

3 K*SMALLK*0MEG*XH*ND0UB*CPCPB*TITLE(8> *M8 *del* alpha, alphaf 
COMMON  /COEFF/  A ( 99) *B (99 ) *C (99 ) .D (99) ,PHI ( 1 0000 ) 

COMMON  /GAMMA/  GAMTE 1 *GAMTE  *PGFF  * GAMFF  *GAMFFS 

COMMON  /STEADY/  PHI S ( 1 0000 ) * AX  IS ( 99) * AX2S ( 99) * BX IS ( 99) * 

1 BX2S(99) *CXS(99) *PHIUBS(99) 

COMMON  /PHI AIR/  CONAFF (500) *PHIAFF (500) 

DIMENSION  PHIOG (99) *0MEGA(99) ,V<99) *EPSGRD(3) 

NAMELIST  /IN/  X * Y * GAMFF *IM*JM*ILE*ITE* JW  *OMEGAH*OMEGAE  *0MEGAP * 

1 EPSGRD*NDUMP*NGRID.NGFF*PGFF,KEPS*NPRINT*SMALLK*IK*XH 
NAMELIST  /CONTRL/  ITAPE 

C 

C TO  START  PROGRAM » STEADY  (PHIS)  DATA  IS  TO  BE  READ  FROM  A DISC 
C FILE  TAPE8.  UNSTEADY  DATA  WILL  NOW  BE  WRITTEN  ON  T APE 7 • 

C 

C TO  RESTART  PROGRAM*  AGAIN  STORE  STEADY  DATA  AS  ABOVE  AND  STORE 
C THE  UNSTEADY  DATA  ON  A DISC  FILE  TAPE9.  NEW  UNSTEADY  DATA  WILL  MOW 
C BE  WRITTEN  ON  TAPE7. 

C 

C READ  STEADY  SOLUTION 

READ  (8)  DUM  *IMS*IM1S*JMS*JM1S*  JWS  »DUM*DUM*ITES*ILES*K  *DEL  * ALPHA  * 
1 DUM*DUM,ND8*DUM,M8*GAM,DUM*DUM*DUM*DUM*DUM*DUM,ALPHAF 
READ  (8)  DUM 
READ  (8)  DUM 

READ  (8)  (DUM* DUM* PH I UBS (I) *I=ILES«ITES) 

READ  (8)  ( AXIS ( I ) * AX2S ( I ) *BX1S(I) *8X2S(I) *CXS(I) *DUM* I=1*IM1S) 
READ  (8)  DUM 
L=IMS* JMS 

READ  (8)  ( PHIS ( I ) ♦ I = 1 *L ) 

C MODIFY  LEADING  AND  TRAILING  EDGE  PHI 
L= ( ILES-1 ) * JMS ♦JWS 
PHIS(L)=.5*(PHIS(L) ♦PHIUBS(ILES) ) 

L=(ITES-1)*JMS*JWS 
PHIS(L)=.5*(PHIS(L) ♦PHIUBS( ITES) ) 

SK=SQRT (ABS(K) ) 

CPCPB=DEL**.6666666667/( ( 1 .♦GAM) *M8**2) «*. 3333333333 
READ  (5*911)  (TITLE(I) *1=1*8) 

READ  ( 5 *CONTRL ) 

IF  (ITAPE. EQ.O)  GO  TO  10 
C READ  DATA  FROM  RESTART  TAPE 

READ  (9)  NITERG*IM*IM1*JM*JM1,JW*JWP1 ♦ JWM1 * ILE* ITE*GAMTE1 * 

1 GAMFF *0MEG*SMALLK*DYBU1*DYBU2*DYBL1*DYBL2*ND0UB*XH 
READ  (9)  (X(I) *I=1*IM) 

READ  (9)  (Y(I) *I=1*JM) 

READ  (9)  (FPU(I) *FPL(I) •PHIUB(I) *I=ILE*ITE) 

READ  (9)  (AX1 (I) , AX2 ( I ) »BX1 (I) *BX2 (I)*CX(I)*DX(I)*I=1*IM1) 

READ  (9)  (AY1 (I) *AY2(I) *DY(I) *1=1* JM1) 

L=IM*JM 
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READ  (9)  (PHI (I) *I=1*L> 

NPT=2* ( IM-3) O*  JM 

READ (9)  (CONAFF(I) *PHIAFF(I> , I=1.NPT> 

GAMTE=GAMTE 1 
gamffs=gamff 

I K = 0 

READ  (5* IN) 

WRITE  (6*900) 

WRITE  (6.901)  NITERG 
NITERG=0 

WRITE  (6*913)  K ♦ CPCP8 

C THE  IK  OPTION  IS  USED  TO  BOOT  STRAP  TO  DIFFERENT  REDUCED  FREQUENCIES 
C AND/OR  MODES  OF  OSCILLATION 
IF  (IK.EQ.O)  GO  TO  15 

0MEG=SMALLK*M8**2/( ( 1 . ♦GAM > *DEL*M8**2 )** .6666666667 
IFAR=0 

IF  ( OMEG.GT . 0 • ) IFAR=1 
CALL  FARFLD  ( IFAR ) 

GO  TO  15 

C START  PROBLEM  FROM  SCRATCH 
10  CONTINUE 
READ  (5. IN) 

N I TERG=0 
NDOUB=0 

0MEG=SMALLK*M8**2/ ( ( 1 . ♦GAM) *DEL*M8«*2) **.6666666667 

GAMTE=GAMFF 

GAMTE 1=GAMFF 

IM1=IM-1 

JM  1 = JM - 1 

JWP1=JW^  1 

JWM 1 = JW- 1 

ERR  = CMPLX ( 0 . ♦ 0 • ) 

DO  3 I = ILE ♦ I TE 
PHIUB ( I ) =ERR 
3 CONTINUE 

IF  (K.LT.O.)  GO  TO  20 

C INITIAL  GUESS  FOR  SUBSONIC  CASE  (INTERIOR  ONLY) 

DO  6 I =3 ♦ IM 1 
M= ( I - 1 ) * JM 
DO  7 J = 2 * JM 1 
L = M ♦ J 

IF  ( Y ( J) .EQ.O . ) GO  TO  8 

PHI (L)=.15915A943l *GAMFF* ( AT  AN ( X ( I)/(SK*Y(J)))*SIGN(1 .570796325* 

1 Y ( J)  ) ) 

GO  TO  7 
8 CONTINUE 

PHI (L)=CMPLX(0.»0.) 

7 CONTINUE 
6 CONTINUE 
GO  TO  30 
20  CONTINUE 

C INITIAL  GUESS  FOR  SUPERSONIC  CASE  (INTERIOR  ONLY) 

DO  21  1=3. IM1 
M=(I-1)*JM 
DO  22  J=2* JM1 
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L=M*J 

SKY=SK*ABS(Y(J> ) 

PHI <L)=CMPLX(0.*0.) 

IF  (X  ( I)  .LT.SKY^XH.OR.X ( I ) .GT.SKY^l . ) GO  TO  22 
PHI (L)=(X(I)-SKY-XH) *S IGN ( 1 . ♦ Y ( J) ) /SK 
22  CONTINUE 
21  CONTINUE 
30  CONTINUE 

c initialize  finite  difference  coefficients  and  farfield 

CALL  INITAL 
IFAR=1 

CALL  FARFLD  ( IFAR) 

15  CONTINUE 

OMEG2I=CMPLX(0.*2.*OMEG) 

WRITE  (6t IN) 

WRITE  (6*900) 

IF  (ITAPE.EQ.O)  WRITE  (6.913)  K.CPCPB 
KGRD=1 

C RE-CYCLE  POINT  FOR  GRID  ITERATION 

50  CONTINUE 

ERROR  = CMPLX ( 0 • « 0 • ) 

NIT=NITERG 
NI TERG=NI TERG* 1 

IF  ( MOD (NI TERG *NPR I NT ) .EQ.O)  CALL  PRINT (NIT) 

IF  (K.LT.O.)  GO  TO  51 

IF  (MOD (NI TERG *NGFF) . NE.O)  GO  TO  51 

GAMFFS=GAMFF 

CALL  GAMFUN 

I F AR  = 0 

CALL  FARFLD  (IFAR) 

WRITE  (6*910)  NITERG.GAMFF.GAMTE 

51  CONTINUE 

DO  1 I =2*  JM 1 
V ( I ) =K 

1 CONTINUE 

C BEGIN  LOOP  ON  GRID 

INCR=2** (NDB-NDOUB) 

I S=3- I NCR 
DO  100  1=3* IM1 
IS=IS* INCR 

C CHECK  FOR  AIR  FOIL 
IFLAG=0 

IF  (ILE.LE.I.AND.I.LE.ITE)  IFLAG=1 

M=(I-1)*JM 

MS= ( IS-1 ) * JMS 

C SAVE  THIS  COLUMN  OF  PHI 
DO  2 J=2* JM1 
L=M  ♦ J 

PHIOG ( J) =PHI (L) 

2 CONTINUE 

C BEGIN  LOOP  ON  COLUMN 
JS=2- I NCR 
DO  150  J=2 ♦ JM 1 
JS= JS* INCR 

C CALCULATE  CELL  INDICES 
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L=M*J 

LR=L*JM 

LL=L-JM 

LLL=LL-JH 

LB=L-1 

LT=L ♦ 1 

C CALCULATE  CELL  INDICES  FOR  PHIO 
LS=MSOS 
LSR=LS* JMS 
LSL=LS-JM5 
LSLL=LSL- JMS 

C CALCULATE  V AND  PHIXX  FROM  STEADY  SOLUTION 

VV=K-AX1S(IS)* (PHIS (LSR)-PHIS(LS) )-AX2S(IS)* (PHIS (LS) -PHIS (LSD ) 

vvs=vv 

IF  (VV.LT.O.)  GO  TO  60 
C ELLIPTIC 

PHI XX=BX1 S( IS)* (PHIS (LSR) -PHIS (LS) > -BX2S < I S > * (PHI S (LS) -PHI S (LSL ) > 
OMEGA (J)=OMEGAE 
GO  TO  199 
60  CONTINUE 

OMEGA (J)=OMEGAP 
IF  (V(J> .GT.O.)  GO  TO  65 
C HYPERBOLIC 

OMEGA (J)=OMEGAH 

VV=K-CXS( IS-1 )*(PHIS(LS)-PHIS(LSL) )-CXS(IS-2)* 

1 (PHIS(LSL)-PHIS(LSLL) ) 

65  CONTINUE 
C PARABOLIC 

PHI XX  = BX IS ( IS-1 >* (PHIS (LS) -PHIS (LSL > >-BX2S( IS-1)* 

1 (PHIS(LSL)-PHIS(LSLL) ) 

199  CONTINUE 
V ( J)=VVS 
TPH I P=PH I (LR) 

TPHIL=PHI (LL) 

IF  ( I .EQ. ILE-1 .AND.J.EQ. JW)  PH  I (LR ) = .5* ( PH  I (L* JM ) ♦PH  I UB ( ILE ) ) 

IF  (I.EQ.ITE*l.ANO.J.EQ.JW)  PH  I ( LL ) = .5* ( PH  I UB ( I TE ) *PH I (L- JM ) ) 

C SET  UP  TRIDIAGONAL  MATRIX  TO  SOLVE  FOR  PHI<I,J) 

C A = PH I ( I » J* 1 ) B = PHI ( 1 1 J)  C = PHI ( I ♦ J-l ) D=CONST  ANT 

IF  < IFLAG.EQ. 1 .AND. J.EQ. JWP1 > GO  TO  250 
IF  ( IFLAG.EQ. I .AND. J.EQ. JWM1 ) GO  TO  270 
IF  ( IFLAG.EQ. 1 .AND. J.EQ. JW)  GO  TO  290 
PART=CMPLX(0.»0.) 

IF  (I.LE.ITE)  GO  TO  201 
C KUTTA  CONDITION 

SIGI=(X(I)-1.)*(GAMFF-GAMTE)/(X(IM1)-1.)*GAMTE 
IF  (J.EQ.JWM1)  PART=.5*AY1 ( J)*SIGI 
IF  (J.EQ.JW)  PART=.5* (AYl ( J)-AY2( J) )*SIGI 
IF  (J.EQ.JWP1)  PART=-.5*AY2( J)*SIGI 
201  CONTINUE 

IF  (VVS.LT.O.)  GO  TO  220 
C ELLIPTIC 

B ( J) =- ( VV» (BX1 ( I ) ♦BX2( I >)*AY1(J)*AY2<J>* (0MEG2I ♦PHIXX)* 

1 ( A X 2 ( I ) -AX  1(1))) 

D ( J ) =- VV* ( 8X 1 ( I ) *PH I <LR> ♦BX2(I)*PHI (LL) ) ♦ ( 0MEG2 1 ♦PH  I XX ) * 

1 ( AX  1 ( I ) *PH I ( LR ) -AX2 ( I ) *PHI (LL) ) ♦PART 
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IF  (J.EQ.2)  GO  TO  202 
IF  (J.EQ.JM1)  GO  TO  203 
A ( J) = AY1 ( J) 

C ( J) =AY2 ( J) 

GO  TO  200 

C BOTTOM  BOUNDARY 

202  CONTINUE 

A ( J) =AY1 ( J) 

D(J)=D<J)-AY2(J)*PHI (LB) 

GO  TO  200 
C TOP  BOUNDARY 

203  CONTINUE 

C ( J) =AY2 ( J) 

D(J)=D( J)-AY1 ( J)*PHI (LT) 

GO  TO  200 

C HYPERBOLIC  AND  PARABOLIC 
220  CONTINUE 

B(J)=VV*BX1 (I-l)-AYl (J)-AY2(J)-( 0MEG2I ♦PH IXX)*CX(I-1) 

D(J)=VV*(BX1 ( I — 1 ) *PHI (LL)+BX2(I-1)* (PHI (LL ) -PH  I (LLL ) ) )- 

1 ( 0MEG2 I ♦PHI XX ) *(CX(I-1)*PHI (LL ) ♦ CX ( 1-2) * (PHI (LL)-PHI (LLL) > > 

2 ♦PART 

IF  (J.EQ.2)  GO  TO  222 
IF  (J.EQ.JM1)  GO  TO  223 
A ( J) =AY1 (J) 

C ( J ) =AY2 ( J) 

GO  TO  200 

C BOTTOM  BOUNDARY 

222  CONTINUE 

A ( J) =AY1 ( J) 

IF  (K.LT.O.)  GO  TO  224 
D(J)=D(J)-AY2(J)*PHI (LB) 

GO  TO  200 

224  CONTINUE 
BT=SK4.5*(K-VV)/SK 

B ( J) =B ( J)-2.*CMPLX  <BT*CX  < I - 1 ) .OMEG/SK ) /DY  < 2 ) 

D(J)=D(J)-2.«BT*(CX(I-1)*PHI (LL ) -CX ( 1-2) * (PHI (LL ) -PH  I (LLL ) ) ) /DY ( 2) 
GO  TO  200 
C TOP  BOUNOARY 

223  CONTINUE 

C ( J) =AY2 ( J) 

IF  (K.LT.O.)  GO  TO  225 
D ( J) =D ( J) -AY  1 ( J ) * P H I (LT) 

GO  TO  200 

225  CONTINUE 
BT=SK^.5*(K-VV)/SK 

B(J)=B(J)-2 • *CMPLX (BT*CX(I-1) .OMEG/SK ) /DY ( J-l ) 
D(J)=D(J)-2.*BT*(CX(I-1)*PHI (LL ) -CX ( 1-2 ) * ( PHI (LL)-PHI (LLL) ) )/ 

1 DY(J-l) 

GO  TO  200 

C BODY  BOUNDARY  ABOVE 
250  CONTINUE 

IF  (VVS.LT.O.)  GO  TO  260 
C ELLIPTIC 

A ( J) =DYBUI 

B(J)=-(VV*(BX1(I).8X2(I))^DYBU1* (0MEG2I ♦PHIXX) * ( AX2 ( 1 ) -AX  1 ( I ) ) ) 
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C ( J > =0 . 

D(J)=DYBU2*FPU( I)-VV*(BX1 (I)*PHI (LR) ♦BX? ( I ) *PHI (LL  > > ♦ (OMEG2I ♦ 

1 PHI  XX ) * (AX  1 ( I ) *PH I (LR)-AX2(I)*PHI(LL) > 

GO  TO  200 

C HYPERBOLIC  AND  PARABOLIC 
260  CONTINUE 
A ( J) =DYBU1 

B ( J) =VV»BX1 <1-1 >-DYBUl-(0MEG2I*PHIXX) *CX ( 1-1 ) 

C ( J) =0. 

D ( J ) =D YBU2*FPU  ( I > ♦ VVMBXl  (I-1)*PHI  (LL > >BX2 ( 1-1 > * (PHI  <LL>- 
1 PHI (LLL ) ) )-(0MEG2I*PHlXX)*(CX(I-l)*PHI (LL ) *CX ( 1-2) * (PHI (LL > - 
1 PHI  (LLL))) 

GO  TO  200 

C BODY  BOUNDARY  BELOW 
270  CONTINUE 

IF  (VVS.LT.0.)  GO  TO  280 
C ELLIPTIC 

A ( J ) =0 . 

B ( J)  =-  (DY8LUW*(RXl  ( I ) *BX2  ( I ) ) ♦ ( 0MEG2 1 ♦PH  I XX ) * ( AX2  ( I ) -AX  1 ( I ) ) ) 
C< J)=DYBL1 

D( J) =-DYBL2*FPL ( I >-VV* (BX1 ( I ) *PHI (LR) 4BX2 ( I ) *PHI (LL) ) ♦ (0MEG2I^ 

1 PH  I XX ) * ( AX  1 ( I ) *PH I ( LR ) -AX2 ( I ) *PH I (LL) ) 

GO  TO  200 

C HYPERBOLIC  AND  PARABOLIC 
280  CONTINUE 
A ( J) =0. 

B(J)=VV*BX1 (I-1)-DYBL1-(0MEG2I^PHIXX)*CX(I-1) 

C ( J ) =D YBL  1 

D (J) =-DYBL2*FPL ( I ) ♦VV* (BX1 (I-1)*PHI (LL)*BX2(I-1)*(PHI (LL) 

1 -PHI  (LLL)  ) ) - (0MPG2I *PHIXX) * (CX ( 1-1 ) *PHI  < LL ) ♦CX  ( I -2 ) * ( PH  I ( LL ) 

2 -PHI (LLL))) 

GO  TO  200 

C BODY  BOUNDARY  J=JW 
290  CONTINUE 
A ( J) =0. 

6 ( J ) =CMPLX ( 1 • ♦ 0 • ) 

C(J)=0. 

D ( J ) =PH I (L) 

200  CONTINUE 

PHI (LR) =TPHIR 
PHI (LL) =TPHIL 
150  CONTINUE 

C TRIDIAGONAL  MATRIX  IS  SET  NOW  SOLVE  FOR  COLUMN  OF  PHI 
CALL  TRI  (I) 

C RELAX  PHI,  FIND  ERROR*  AND  MOVE  TO  NEXT  COLUMN 
DO  295  J = 2 * JM 1 
L = M*  J 

ERR=OMEGA(J)*(PHI (L)-PHIOG(J) ) 

PHI (L)=PHIOG(J) ♦ERR 

IF  (CABS (ERR) *LT.CABS (ERROR) ) GO  TO  295 
ERROR=ERR 
LERROR=L 
295  CONTINUE 

IF  ( IFLAG.NE . 1 ) GO  TO  100 
L=M ♦ JW 
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c 


c 


PHI <L)=PHI (L-l ) ♦DY  ( JWMl ) * (PHI (L-l)-PHI (L-2) >/0Y ( JW-2) 
PHIUB ( I ) =PHI (L* 1 ) -DY ( JW)  * (PH  I <L*2)-PHI < L ♦ 1 ) ) /OY ( JWP1 ) 
IF  (I.EQ.ITE)  6AMTE=PH IUB ( I ) -PHI (L ) 

IF  (CABS(GAMFF) .LE.l .E-8)  GAMTE=CMPLX ( 0 . * 0 . ) 

IF  (K.LT.O.)  GAMFF=GAMTE 

IF  (I.EQ.ITE. AND. NITERG.EQ.I)  GAMTE1=GAMTE 
100  CONTINUE 

PRINT  OUT  ERROR  AFTER  EACH  GRID  SWEEP 
WRITE  (6*905)  N I TERG* ERROR  *L ERROR 
IDOU8=0 

IF  (CABS (ERROR) .LE .EPSGRD (KGRD) ) GO  TO  300 
IF  (NITERG.EQ.NGRID)  GO  TO  310 
IF  (MOD(NITERG*NDuMP) .EQ.0)  GO  TO  310 
GO  TO  50 

300  CONTINUE 
KGRD=KGRD* 1 
IDOUB=l 

GO  TO  310 

301  CONTINUE 

IF  (K.LT.O.) 

GAMFFS=GAMFF 
CALL  GAMFUN 
WRITE  (6*910) 

302  CONTINUE 
CALL  FPRINT 
WRITE  (6*900) 

WRITE  (6*906) 

CALL  DOUBLE 
WRITE  (6*914) 

(6*902) 

(6*903) 

(6*904) 

(6*903) 

50 


GO  TO  302 


niterg*gamff*gamte 


niterg 


( X ( I ) *I  = 1*IM> 
(Y(I) *I=1*JM) 


IM*JM*JW*ILE*ITE 

WRITE 
WRITE 
WRITE 
WRITE 
GO  TO 
310  CONTINUE 
TAPE  DUMP 

WRITE  (7)  NITERG,IM,IM1,JM*JM1*JW*JWP1*JWM1,ILE*ITE.GAMTE1* 
1 GAMFF*0MEG*SMALLK*DYBU1 *DYBU2*DYBL1*DYBL2*ND0UB*XH 
write  (7)  ( X C I ) *I  = 1*IM) 

(Y(I)*I=1*JM) 

(FPU (I), FPL (I) ♦ PH IUB (I) *I=ILE*ITE) 

(AX1 ( I ) , AX2  < I ) *BX1 (I)*BX2(I) *CX(I)*DX(I) , I = 1 ♦ I M 1 ) 
(AY1 (I) ,AY2(I) *DY(I) , I=1*JM1) 


WRITE 
WRITE 
WRITE 
WRITE 
L=IM*JM 
WRITE  (7) 


(7) 

(7) 

(7) 

(7) 


(PHI  < I ) * 1 = 1 *L) 


NPT  = 2* ( IM-3) >3*JM 


WRITE  (7)  (CONAFF ( I ) *PHI AFF (I) *1=1* NPT) 
END  FILE  7 

WRITE  (6*907)  NITERG 
CALL  PRINT (NITERG) 

IF  (KGRD.GT.KEPS)  GO  TO  320 
IF  (NITERG.EQ.NGRID)  GO  TO  330 
IF  ( IDOUB.EQ. 1 ) GO  TO  301 
GO  TO  50 
320  CONTINUE 
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WRITE  (6*908)  NITERG 
GO  TO  350 
J30  CONTINUE 

WRITE  (6*909)  NITERG 

900  FORMAT  (1H1) 

901  FORMAT  (IH  */**  CASE  IS  BEING  RESTARTED  AT  I TERAT I ON* 15 ) 

902  FORMAT  (1H  */♦*  X ( I > * 1 = 1 ♦ IM* ) 

903  FORMAT  (10E13.5) 

904  FORMAT  (1H  ./**  Y ( I ) » 1=1 ♦ JM*) 

905  FORMAT  (1H  ♦ /*«  AT  ITERATI0N*I5*  THE  MAXIMUM  ERROR  =*2E13.5 

1 * AND  OCCURRED  AT  N0DE*I5> 

906  FORMAT  (lH  ♦ /♦*  THE  NUMBER  OF  NODES  IS  BEING  DOUBLED  AT  ITERATION* 

1 15) 

907  FORMAT  <1H  */,*  TAPE  HAS  BEEN  DUMPED  AT  ITERAT ION* 1 5) 

908  FORMAT  (1H  */**  SOLUTION  HAS  CONVERGED  TO  DESIRED  ACCURACY  AT  ITER 
1 AT I0N*I5) 

909  FORMAT  C 1 H ./♦*  MAXIMUM  NUMBER  OF  ITERATIONS  HAS  BEEN  REACHED*  CAS 
IE  IS  BEING  TERMINATED  AT  I TERAT I ON* 1 5) 

910  FORMAT  (1H  ♦ /♦*  UPDATE  GAMFF  AND  FARF1ELD  AT  ITERATI0N*I5 

1 * GAMFF  =*2E13.5*  GAMTE  =*2E13.5) 

911  FORMAT  (8A10) 

913  FORMAT  (1H  ♦ /,*  SIMILARITY  PARAMETER  (K)  =*E13.5*/**  SCALING  FACTO 
1R  ( CP/CPBAR ) =*E13.5) 

914  FORMAT  (lH  ♦/**  IM  = *I4*  JM  = *I4*  JW  = *I4*  ILE  =*I4*  ITE  =*I4) 

350  CONTINUE 

CALL  FPRINT 
END 

SUBROUTINE  DOUBLE 

COMPLEX  PHIUB*B*D*PHI*GAMTE1 * GAMTE  * GAMFF  * GAMFFS ♦ S I G I »PDUM 
REAL  K ♦ M8 

COMMON  /DELTA/  DX ( 99) *DY ( 99) * AX  1 (99) * AX2 (99) »BX 1 ( 99 ) *BX2 (99 ) , 

1 CX ( 99 ) *AY1 (99) *AY2<99) *X(100) » Y ( 1 00 > *FPU < 99 > *FPL(99) *PHIU6(99) ♦ 

2 IM*IM1*JM»JM1*JW* JWP 1 ♦ JWM 1*ILE*ITE*DYBU1*  DYBU2*DYBL 1 ♦ DYBL2 » 

3 K»SMALLK.0MEG*XH»ND0UB*CPCPB*TITLE(8) *M8 ♦ DEL  * ALPHA » ALPHAF 
COMMON  /COE FF/  A (99) »B (99) ,C (99) *D (99) *PHI ( 10000) 

COMMON  /GAMMA/  GAmTE 1 » GAMTE  tPGFF ♦ GAMFF ♦ GAMFF S 
ND0UB=ND0UB>1 
C CHANGE  GRID  INDICES 

I MAX  = 2* I M— 4 
IMAXl=IMAX-l 
JMAX=2* JM-3 
JM AX  1 = JMAX- 1 
JWN=2* JW-2 
JWNP 1 = JWN* 1 
JWNM 1 = JWN- 1 

I TEN=2* I TE-3 

DIST=l./(l.*DY(JWPl)/DY(JW) ) 

IK=ITE 
C CHANGE  X 
L=  IM 

DO  30  I =5 ♦ I MAX  1 *2 

I I = IMAX 1-1 ♦S 
L=L- 1 

X < 1 1 ) =X (L ) 

X < 1 1 — 1 )=.5*(X(L)*X ( L — 1 ) ) 
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30  CONTINUE 

X(IMAX)=2.*X(IMAX1)-X(IMAX1-1) 

OELT=X (4) -X ( 3) 

X ( 2 ) =X ( 3) -OELT 
X ( I ) =X (2) -OELT 
C FIND  NEW  LEADING  EDGE 
DO  31  1=3  ♦ I MAXI 
IF  (X(I) .GE.O.)  GO  TO  32 

31  CONTINUE 

32  ILEN= I 
C CHANGE  Y 

L=JM 

DO  40  I=2» JHAXl *2 
1 1= JMAX1-I *2 
L=L- 1 

Y < 1 1 ) =Y (L ) 

Y (II-1)  = #5*<Y (L) ♦ Y < L — 1 ) ) 

IF  (II-1.EQ.JWNP1)  Y < II-l ) =Y  (L-1)*DIST*(Y(L)-Y(L-1) ) 

IF  (II-1.EQ.JWNM1)  Y(II-1)=Y(L>*DIST*(Y(L-1>-Y(L) > 

40  CONTINUE 

Y< JMAX)=2.*Y( JHAX1)-Y( JMAX1-1) 

Y(1)=2.*Y(2)-Y(3) 

C MOVE  THE  PHI(S)  < INTERIOR  ONLY)  AND  INTERPOLATE  FOR  PHI  (I=ODD> 
L=IM1* JM 

DO  10  I=3»IMAXlt2 
I I = IMAXl-I *3 
LL= ( 1 1 -1 ) * JMAX 
DO  15  J=2t JMAXl *2 
JJ=JMAX1-J*2 
LLL=LL*JJ 
L=L- 1 

PHI (LLL ) =PHI (L) 

PHI (LLL-1 ) =.5* (PHI (L) ♦PHI (L-l) > 

IF  (JJ.NE.JWN)  GO  TO  15 

PHI (LLL-1 ) =PHI (L) *DIST* (PHI (L-l) -PHI (L) ) 

IF  (II.GT.ITEN)  GO  TO  16 
IF  (II.GE.ILEN.AND.II.LE.ITEN)  GO  TO  17 
PHI (LLL* 1 ) =PHI (L) *DIST* (PHI (L*1)-PHI (L) ) 

GO  TO  15 

16  CONTINUE 

SIGI= (X  < 1 1 ) -1 . ) * (GAMFF-GAMTE) / (X ( IMAX1 ) -1 . ) *GAMTE 
PDUM=PHI (L) *.5*SIGI 

PHI (LLL* 1 ) =PDUM*DIST* (PHI (L*1)-PDUM) 

PDUM=PHI (L)-.5*SIGI 

PHI (LLL-1 ) =PDUM*DIST* (PHI (L-l)-PDUM) 

GO  TO  15 

17  CONTINUE 

PHI (LLL* 1 ) =PHIUB (IK)*DIST*(PHI (L* 1 ) -PHIUB ( IK ) ) 

IK=IK-1 
15  CONTINUE 
L=L-2 

10  CONTINUE 

C FILL  IN  OTHER  PHI (S)  (I=EVEN) 

IMAX2=IMAX 1-1 
DO  20  I=4»IMAX2*2 
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I  I = IMAX2-I*4 
L= ( 1 1 - 1 ) * JMAX 
DO  25  J=2  , JMAX 1 
JJ=JMAX1-J*2 
LL=L* J J 

PHI (LL) =.5* (PHI <LL* JMAX) ♦PHI (LL-JMAX) ) 

25  CONTINUE 
20  CONTINUE 
IM=IMAX 
IM1=IMAX1 
JM= JM AX 
JM 1 = JMAX 1 
JW= JWN 
JWP 1 = JWNP 1 
JWM 1 = JWNM 1 
ITE=ITEN 
ILE=ILEN 

KTE=(ITE-l)*JMOW 

PHI (KTE* JM) =.5*PHI (KTE) ♦.25»GAMTE* .5*PHI <KTE»2*JM) 

C INITIALIZE  FINITE  DIFFERENCE  COEFFICIENTS  AND  FARFIELD 
CALL  INITAL 
L=( ILE-1)*JM>JW 
DO  50  I=ILE,ITE 

PHIUB(I)=PHI <L*1)-DY<JW>*<PHI <L^2)-PHI (L^l) >/DY(JWPl> 

L=L ♦ JM 
50  CONTINUE 

L=(ILE-2)*JM»JW 

PHI  (L)  =.5*  (PHI  (L-JM) ♦.5*<PHIU6(ILE) ♦PHI  (LOM)  ) ) 

I F AR  = 1 

CALL  FARFLD  ( IFAR ) 

RETURN 

END 

SUBROUTINE  FARFLD  (IFAR) 

COMPLEX  PHIUB,B,D,PHI ,GAMTE1 , GAMTE , GAMFF ♦GAMFFS , PART  1 , PART2 ♦ 

1  PART3,PART4,PART5,PART3N»H1 , H2» CONST* CONAFF , PHI AFF 
REAL  K»KH»M8 

COMMON  /DELTA/  DX ( 99) , DY ( 99 ) ♦ AX  1 (99 > » AX2 (99) »BX 1 < 99 ) *BX2 ( 99) ♦ 

1 CX(99) »AY1 (99) »AY2(99>  * X < 1 00 > * Y ( 1 00 ) ,FPU < 99 ) ,FPL (99) ♦PHIUB (99) ♦ 

2 IM» IM1 * JM, JM1 * JW» JWP1 ♦ JWM1 , ILE»ITE»DYBU1 »DYBU2  »DYBL 1 » DYBL2  » 

3 K.SMALLK»0MEG*XH»ND0UB»CPCP8» TITLE (8) ,M8,DEL, ALPHA, ALPHAF 
COMMON  /COEFF/  A ( 99 > *B ( 99) ,C ( 99 > * D ( 99 ) ,PH I ( 1 0000 > 

COMMON  /GAMMA/  GAMTE 1 , GAMTE , PGFF »GAMFF ♦ GAMFF S 
COMMON  /PHIAIR/  CONAFF (500) »PHI AFF (500) 

DIMENSION  EX(16)*WT(16) 

DATA  (EX  (I) * 1 = 1 ♦ 16)/. 0483076657 ». 1444719616, .2392873623, 

1 .33 18686023*. 421 351 2761 ,.5068999089,. 5877157572, .6630442669, 

2 .7321821 187, .794483796 *.8493676 137, .896321 1558, .9349060759, 

3 .9647622556, .9856115115, .9972638618/ 

DATA  (WT ( I ) , 1=1 , 16) /. 0965400885, .0956387201 , .0938443991 ♦ 

1 .091 1738786,. 087652093,. 08331 19242,. 0781 938958,. 0723457941, 

2 .0658222228,.0586840935,.0509980593, . 042835898 ,. 0342738629 , 

3 .0253920653, .01 62743947 ,. 0070 1861/ 

DATA  RSTAR  /10./ 

IF  (K.LT.O.)  GO  TO  90 
C SUBSONIC  FARFIELD 
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IF  (IFAR.EQ.O)  60  TO  60 
SK=SQRT(K) 

KH=OMEG/K 

IF  (KH.GE..05)  GO  TO  30 
LFF =0 

DO  10  1=1 » IM 
H= ( I — 1 » * JM 
J3= JM1 

IF  (I.LE.2.0R.I.EQ.IM)  J3=l 
DO  15  J=1 « JMt J3 
L=M*J 
LFF  =LFF ♦ 1 

CONAFF (LFF ) =CMPLX ( 0 • * 0 • ) 

PHI AFF (LFF ) =CONAFF (LFF) 

IF  (Y(J).EQ.O.)  GO  TO  20 

PHI (L) =. 1591549431 *GAMFF« (AT AN (X  < I)/(SK*Y(J)))*SIGN( 1.570796325* 
1 Y ( J ) ) ) 

GO  TO  15 
20  CONTINUE 

PHI (L ) =CMPLX ( 0 • * 0 . ) 

15  CONTINUE 
10  CONTINUE 
RETURN 
30  CONTINUE 

CONST  = *25*0MEG*GAMFF /SK 

LFF  = 0 

ZIP=0. 

PART2=CHPLX ( 0 • » 0 * ) 

PART4=PART2 
DO  35  I = ILE ♦ I TE 

PART1= (PHIUB ( I )-PHI ( (I-l)*JMOW) ) *CEXP (CMPLX (0. *-KH*X ( I ) ) ) 
PART4=PART4*,5* (PARTI 4PART2) * (X ( I ) -ZIP) 

ZIP=X(I) 

PART2=PART  1 
35  CONTINUE 

DO  40  I = 1 * I M 
M=(I-1)*JM 
J3= JM1 

IF  (I.LE.2.0R.I.EO.IM)  J3=l 
DO  45  J=1 ♦ JH* J3 
LFF=LFF ♦ 1 
L=M»  J 

IF  (Y(J).EQ.O.)  GO  TO  46 
PART3=CMPLX ( 0 • • 0. ) 

R=KH*SORT ( (X ( I )-l .) **2*K*Y ( J)**2) 

CALL  HANKEL  (R»H1,H2) 

CONAFF ( LFF ) =H 1 *CMPLX ( 0. ♦ -.25*SK*KH**2*Y(J)/R)*CEXP( CMPLX (0.* 

1 KH*X ( I ) ) ) 

PHI AFF (LFF ) =PART4*C0NAFF (LFF ) 

PARTI =-CEXP ( CMPLX (0.*KH*(X(I)-1.) ) ) *H1/R 
RI=RSTAR 

IF  (R.GT.RSTAR)  RI=R 
S=ABS(KH*(X(I)-1.) ) >RI 
SS=SQRT (S) 

G=S**2*K*(KH*Y(J) )**2 
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GM=S**2-K* (KH*Y< J) ) **2 
AA=GM/(SS*G*»1.5> 

BB= ( 1 .5*SS*.5*K*(KH*Y ( J) >**2/(SS*S> ) /G**l .5-3.*SS*GM/G*»2.5 
CC=3.75*SS*GM/G**2.5 

DD=(9.375*S*SS-1.875*K*(KH*Y( J) ) **2/SS) /G**2.5-18. 75*S*SS* 

1 GM/G**3 . 5 

PART2=1 . 128379 167*CEXP (CMPLX ( 0. . ,7853981635-S) ) *CHPLX (BB*CC» AA-DD) 
PART5=CHPLX(0.*0.) 

IF  (R.GE.RSTAR)  GO  TO  47 
R I =R 

RMID=.5*(RSTAR*RI) 

RAD=.5*(RSTAR-RI) 

SN=1. 

C GAUSSIAN  QUADRATURE 
DO  41  I J=1 *32 
I A=IABS ( I J-16) 

IF  (IJ.LE.16)  IA=IA*1 
IF  (IJ.GT.16)  SN=- 1 • 

RI=RMID»RAD*EX(IA)*SN 
CALL  HANKEL  (RI.H1.H2) 

PART3N=CEXP (CMPLX (0. » -SORT ( R I **2-K* (KH*Y ( J> )**2> ) ) *H2/RI 
PART5=PART5*WT ( I A) *PART3N 

41  CONTINUE 
PART5=PART5*RAD 

47  CONTINUE 

IF  < X ( I ) .LE.l.)  GO  TO  44 
RI=KH*SK*ABS(Y( J) ) 

RM I D= . 5* ( R>R I ) 

RAD=.5*(R-RI) 

SN=  1 • 

C GAUSSIAN  QUADRATURE 
DO  42  I J= 1 ♦ 32 
IA=IABS< IJ-16) 

IF  (IJ.LE.16)  I A= I A ♦ 1 
IF  (IJ.GT.16)  SN  = - 1 . 

RI=RMID>RAD*EX ( I A) *SN 
CALL  HANKEL  (RI,H1*H2) 

PART3N=SIN (SORT (RI**2-K* (KH*Y ( J) >«*2> )*H2/RI 
PART3=PART3*WT ( I A ) *PART3N 

42  CONTINUE 

PART3=PART3*CMPLX (0. *-2.*PAD) 

44  CONTINUE 

PHI (L) = (PARTI ♦PART2*PART3 ♦PARTS) *CONST*Y ( J) ♦PHIAFF (LFF) 

GO  TO  45 
46  CONTINUE 

PHI ( L ) =CMPLX ( 0 • * 0 . ) 

CONAFF (LFF)=PHI (L) 

PHIAFF (LFF ) =PHI (L) 

45  CONTINUE 
40  CONTINUE 

RETURN 
60  CONTINUE 

CONST=GAMFF /GAMFFS 

LFF  = 0 

ZIP=0. 
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KH=OMEG/K 

PART2=CMPLX(0.,0.) 

PART3=PART2 
DO  80  I=ILE» ITE 

PART1= (PHIUB ( I ) -PHI ( ( 1-1 ) *JM* JW) ) *CEXP <CMPLX (0, *-KH*X ( I ) ) ) 
PART3=PART3*.5*(PART1*PART2)*(X(I)-ZIP) 

ZIP=X(I) 

PART2=PART 1 
80  CONTINUE 

DO  70  1 = 1 * IM 

H=(I-1)*JM 

J3=JM1 

IF  (I.LE.2.0R.I.EQ.IM)  J3=l 
DO  75  J=1,JM,J3 
L=M*  J 
LFF=LFF ♦ 1 

PHI (L)=PHI (L)-PHIAFF (LFF) 

PHI AFF (LFF) =PART3*C0NAFF (LFF) 

PHI (L)=PHI (L)*CONST»PHIAFF(LFF) 

75  CONTINUE 
70  CONTINUE 
RETURN 

C SUPERSONIC  FARFIELD 

90  CONTINUE 

DO  91  1=1, IM 

H=(I-1)*JM 

J3=JM1 

IF  (I.LE.2.0R.I.EQ.IM)  J3=l 
DO  92  J=1 , JM, J3 
L=M*  J 

PHI (L ) =CMPLX ( 0 • * 0 • ) 

92  CONTINUE 

91  CONTINUE 
RETURN 
END 

SUBROUTINE  FPRINT 

COMPLEX  phiub,b*d,phi,gamtei,gamte,gamff,gamffs,clift,cmom* 

1  CHINGEtCl *C2*C3»C10*C20»C30*CIRLIF  *PART 
REAL  K»M8 

COMMON  /DELTA/  DX (99) ,DY (99) , AX1 (99) , AX2 (99) ,BX1 (99) »BX2 (99) , 

1 CX (99) , AY1 (99) ,AY2 (99) ,X ( 100) ,Y ( 100) ,FPU(99) , FPL (99) ,PHIUB(99) * 

2 IM » I M 1 , JM » JM 1 , JK  » JWP 1 , JWM1 »ILE,ITE,DYBUl ,DYBU2  »DYBL 1 »DYBL2  * 

3 K,SMALLK,0MEG,XH,ND0UB,CPCPB,TITLE(8) ,M8, del, alpha, ALPHAF 
COMMON  /COEFF/  A (99) ,B (99) ,C (99) ,D (99) ,PHI ( 10000) 

COMMON  /GAMMA/  GAMTE1 ,GAMTE«PGFF,GAMFF,GAMFFS 
CPDEL=CPCPB/DEL 

PART=.5*( (X(ITE>1)-1.)*(GAMFF-GAMTE)/(X(IM1)-1.)*GAMTE) 

I JW=ITE*JM» JW 
PHI ( I JW) =PHI ( I JW ) -PART 
L= ( ILE-2) * JM* JW 
PHIUB ( ILE-1 ) =PHI (L) 

PHIUB ( ITE*1 ) =PHI ( I JW) *2,*PART 
C COMPUTE  CP  LOWER  (B)  AND  CP  UPPER  (D) 

DO  10  I = ILE»  ITE 
M=(I-1)*JM 
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c 


L=M> JW 

B(I)s-2,*(AXl (I)«(PHI (L4 JM) -PHI (L ) ) ♦ AX2 ( I ) • ( PH  I <L ) -PHI (L-JM) ) ) 
1 *CPOEL 

D ( I>=-2.* (AX1 ( I)*(PHIUB ( I ♦ 1 ) -PHIUB ( I ) ) ♦AX2(I) * <PHIUB< I)- 
1 PHIUB(I-l) ) )*CPDEL 
10  CONTINUE 

PHI ( I JW) = PHI ( I JW) ♦PART 
COMPUTE  UNSTEADY  FORCE  COEFFICIENTS 
CIO =B (ILE) -D (ILE) 

CL IFT=C 1 0*X ( ILE ) 

C20=CLIFT 

CM0M=.5*CLIFT*X(ILE) 

C30=CMPLX (0. *0. > 

CHINGE=C30 

IF  (XH.GE.X(ILE) ) GO  TO  25 
C30=C10*(X(ILE)-XH) 

CHINGE=.5*C30*(X(ILE)-XH) 

25  CONTINUE 
I LE 1 = ILE ♦ 1 
DO  30  I = I L E 1 ♦ ITE 
C 1 =B ( I ) -D  < I ) 

C2=C 1 *X ( I ) 


30 


40 


IF  (X(I).GT.XH)  C3=C1MX<I)-XH) 
CLIFT=CLIFT^.5*(C1>C10)*DX(I-1) 
CMOM=CMOM4 .5* (C2^C20) *DX ( 1-1 ) 

DXX  = DX  <1-1 ) 

IF  (X ( I ) .GT.XH.AND.X ( 1-1 ) .LE.XH) 

CHINGE=CHINGE^.5* (C3^C30)*DXX 

C10=C1 

C20=C2 

C30=C3 

CONTINUE 

CIRLIF=2.*GAMTE*CPDEL 
WRITE  (6,900) 

WRITE  (6,901)  (TITLE(I) ,1=1,8) 

WRITE  (6,902)  M8 

WRITE  (6,903)  K 

WRITE  (6,904)  DEL 

WRITE  (6,922)  ALPHA 

WRITE  (6,905)  SMALLK 

WRITE  (6,906)  OMEG 

WRITE  (6,907)  XH 

WRITE  (6,923)  ALPHAF 

WRITE  (6,908)  CPCPB 

WRITE  (6,909)  GAMTE 

GO  TO  40 
CIRLIF 


IF  (XH.GT.O.) 
WRITE  (6,910) 
WRITE  (6,911) 
WRITE  (6,913) 
WRITE  (6,914) 
GO  TO  45 
CONTINUE 
WRITE  (6,921) 
WRITE  (6,912) 
WRITE  (6,913) 


CLIFT, CMOM,CHINGE 


CIRLIF 


CLIFT, CM0M,CHINGE 


DXX=X ( I )-XH 
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WRITE  (6*915) 

45  CONTINUE 

WRITE  (6*916) 

WRITE  (6*917)  < X ( I ) *I  = ILE*ITE) 

WRITE  (6*918) 

WRITE  (6*919)  (D(I) *I=ILE,ITE) 

WRITE  (6*920) 

WRITE  (6*919)  (B  < I ) » I = ILE* ITE) 

900  FORMAT  (1H1) 

901  FORMAT  (30X*8A10) 

902  FORMAT  ( 1 H */*lH  ,/,lH  ♦/,*  MACH  NUMBER  =*E13.5> 

903  FORMAT  (*  SIMILARITY  PARAMETER  (K)  =*E13.5> 

904  FORMAT  (*  THICKNESS  RATIO  =*E13.5) 

905  FORMAT  (*  REDUCED  FREQUENCY  (BASED  ON  CHORD)  =*E13.5> 

906  FORMAT  (*  SCALED  FREQUENCY  (OMEGA)  =*E13.5) 

907  FORMAT  <*  HINGE  POINT  =*E13.5) 

908  FORMAT  (*  CP  SCALING  FACTOR  (CP/CPBAR)  =*E13.5> 

909  FORMAT  (*  SCALED  AIRFOIL  CIRCULATION  =*2E13.5> 

910  FORMAT  (*  LIFT  COEFFICIENT  BASED  ON  CIRCULATION  (PER  UNIT  PITCH  AN 
1GLE  IN  RADIANS)  =*2E13.5) 

911  FORMAT  (1H  ./*1H  ♦ /**  UNSTEADY  FORCE  COEFFICIENTS  (PER  UNIT  PITCH 
1ANGLE  IN  RADIANS)*) 

912  FORMAT  (1H  ♦/♦ 1H  ,/♦*  UNSTEADY  FORCE  COEFFICIENTS  (PER  UNIT  FLAP  A 
INGLE  IN  RADIANS)*) 

913  FORMAT  (1H  */*3X*LIFT  =*2E13.5*/*3X*M0MENT  ABOUT  (X=0>  =*2E13.5*/* 
1 3X*HINGE  MOMENT  =»2E13.5) 

914  FORMAT  (1H  */*lH  ,/♦*  PRESSURE  COEFFICIENTS  ON  THE  AIRFOIL.  UNSCAL 
1ED  (PER  UNIT  PITCH  ANGLE  IN  RADIANS)*) 

915  FORMAT  (1H  */*lH  */♦*  PRESSURE  COEFFICIENTS  ON  THE  AIRFOIL.  UNSCAL 
1ED  (PER  UNIT  FLAP  ANGLE  IN  RADIANS)*) 

916  FORMAT  (1H  */*3X*AIRF0IL  COORDINATE  =*) 

917  FORMAT  (3XE13.5* 13XE13.5* 13XE13.5* 13XE13.5* 13XE13.5) 

918  FORMAT  (1H  */*3X*AIRFOIL  PRESSURE  COEFFICIENTS*  UPPER  =*> 

919  FORMAT  (3X10E13.5) 

920  FORMAT  (1H  ./ * 3X*A IRFO IL  PRESSURE  COEFFICIENTS*  LOWER  =*) 

921  FORMAT  (*  LIFT  COEFFICIENT  BASED  ON  CIRCULATION  (PER  UNIT  FLAP  ANG 
1LE  IN  RADIANS)  =*2E13.5) 

922  FORMAT  (*  MEAN  AIRFOIL  ANGLE  OF  ATTACK  (RADIANS)  =*E13.5> 

923  FORMAT  (*  MEAN  FLAP  ANGLE  (RADIANS)  =*E13.5) 

RETURN 

END 

SUBROUTINE  GAMFUN 

COMPLEX  GAMTE1*GAMTE*GAMFF*GAMFFS 

COMMON  /GAMMA/  GAMTE1 *GAMTE*PGFF.GAMFF*GAMFFS 

GAMFF =GAMTE 1 *PGFF* (GAMTE-GAMTE 1 > 

GAMTE 1 =GAMTE 

RETURN 

END 

SUBROUTINE  HANKEL  (R.H1.H2) 

COMPLEX  HO *H1 * H2 

IF  (R.GT.3.)  GO  TO  10 

X2=(R/3.)**2 

X4=X2*X2 

X6=X4*X2 

X8=X6*X2 
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X 1 0=X8*X2 
X 12=X 1 0*X2 

A0=1 .-2. 2499997*X2*1 ,2656208*X4-. 31 63866*X64.0444479»X8-. 0039444 
1 *Xl0*.00021*Xl2 

B0=2./3. 14 159265* ALOG (.5*R)*A0*.  36746691 ♦ • 6 0559366* X 2- . 74350384 
1 *X4*.253001 17*X6-. 0426 12 14*X8*. 0042791 6*X10-.00024846*X 12 

A1=R*(.5-.56249985*X2*.21093573*X4-.03954289*X6».00443319*X8 
1 -.00031761*X10*.00001109*X12> 

Bl=(2./3. 14159265*R*AL0G ( ,5*R) *A1-.6366198* .22 12091 *X2»2. 1682709 
1 *X4-1 ,3164827*X6*.3123951*X8-.0400976*X104.0027873*X12)/R 

GO  TO  20 
10  CONTINUE 
XI 1=3. /R 
X 1 2 = X 1 1 *X  1 1 
X 1 3 = X 1 2*X 1 1 
X I4  = X I 3*X 1 1 
X 1 5 = X 1 4*X 1 1 
XI6  = XI5*XI  1 

F =.79788456-. 0000 00 77*X I 1-.0055274*XI2-.00009512*X 1 3* .00137237 
1 *Xl 4- .0007280 5* xIS^.0001 4476* XI6 
TH=R-. 785 398 16- .04 166397* XI 1-.00003954*XI2*.00262573*XI3 
1 -. 000541 25*X I4-. 00029333*X 15*. 000 13558*X 16 

AO=F*COS ( TH) /SORT (R) 

BO=F*SIN(TH)/SQRT(R) 

FI* • 7 9788456 ♦ ,00000156*XI1*.01659667*XI2*.00017105*XI3-.00249511 
1 *Xl4^.0011 3653* X I5-.00020033*XI6 

THl=R-2. 356 19449 ♦ . 124996 1 2*X 1 1 ♦ • 0000565*X 12- .00637879*XI3 
1 ♦.00074348*XI4*.00079824*XI5-.00029166*XI6 

A1=F1*C0S(TH1)/SQRT(R) 

B1=F1*SIN(TH1 )/SQRT (R) 

20  CONTINUE 

H0=CMPLX (A0*-B0> 

H1=CMPLX (A1 *-Bl ) 

H2=2.* H1/R-H0 

RETURN 

END 

SUBROUTINE  INITAL 
COMPLEX  PHIUB 
REAL  K ♦ MB 

COMMON  /DELTA/  DX (99) *DY (99) » AX1 (99) » AX2 (99) *BX1 (99) »BX2 (99) * 

1 CX(99) »AY1 (99) ♦ A Y2 ( 99 ) » X ( 1 00 ) • Y ( 1 00 ) * FPU (99) * FPL (99) ♦ PHIUB (99) ♦ 

2 IM» IM1 ♦ JM* JM1 * Jw* JWP1 ♦ JWM1 » I LE ♦ I TE » DYBU 1 ♦ D YBU2 * DYBL 1 *DYBL2 ♦ 

3 K.SMALLK.0MEG*XH*ND0UB«CPCPB,TITLE(8) ,M8»DEL»ALPHA»ALPHAF 
C CALCULATE  DX 

DO  15  1=1 ♦ IM1 
DX ( I ) =X ( I ♦ 1 ) -X ( I ) 

15  CONTINUE 
C CALCULATE  DY 

DO  25  1 = 1*  JM 1 
DY  ( I ) =Y ( I ♦ 1 ) -Y ( I ) 

25  CONTINUE 

DO  30  1=2* IM1 

AX1(I)=0X(I-1)/(DX(I)*(DX(I-1)*DX(I) ) ) 

AX2(I)=DX(I)/(DX(I— 1)*(DX(I— 1)*DX(I) ) ) 

BX1  ( I ) =2 . *AX 1 ( I )/DX ( 1-1 ) 
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BX2(I)=2.*AX2(I)/DX(I) 

CX ( I ) =#5/DX ( I ) 

30  CONTINUE 

CX ( 1 ) = #5/DX ( 1 ) 

DO  40  1=2 » JM1 

AY1  (I)=2./(DYU)*(DY(I)  ♦DY(I-l) ) ) 
AY2(I)=2./(DY(I-I)*(DY(I>*DY(I-1) )> 

40  CONTINUE 

IF  (K.GT.O.)  GO  TO  50 

C TOP  AND  BOTTOM  ANO  RHS  BOUNDARY  CONDITIONS 
AX1 (IM1)=0. 

AX2 ( IM 1 ) =0 • 

BX1 (IM1)=0. 

BX2(IM1)=2./DX(IM1-1)**2 
AY1 (2)=2./DY<2)**2 
AY2(2)=0. 

AY1 (JMI)=0. 

AY2( JM1)=2./DY(JM1-1)**2 
50  CONTINUE 

DYBU1=2»/ ( (DY ( JWP]  > ♦2.*DY ( JW) ) *DY ( JWP1 ) ) 

DYBU2=DY  ( JWP1 > *DYBU1 

DYBLl=2./( (DY < JW-2) ♦2,*DY ( JWM1 ) ) *DY ( JW-2) ) 

DYBL2=DY< JW-2)*DYBLl 
C SET  AIRFOIL  BOUNDARY  CONDITION 
DO  45  I=ILE* ITE 
FPU ( I ) =0 • 

IF  (X(I).GE.XH)  FPU ( I ) =- 1 • 

FPL(I)=FPU(I) 

45  CONTINUE 
RETURN 
END 

SUBROUTINE  PRINT  (NITERG) 

COMPLEX  PHIUB  »B  *D*PHI *GAMTE 1 *GAMTE .GAMFF  *GAMFFS*PART 
REAL  K *M8 

COMMON  /DELTA/  DX (99) *DY (99) ♦ AX1 (99) * AX2 (99) *BX1 (99) *BX2 (99) * 

1 CX (99) * AY  1 (99) »AY2(99) *X  ( 100) *Y ( 100) *FPU(99) ♦FPL(99) *PHIUB(99) » 

2 IMvIMl* JMt JM1*JW*JWP1* JWM1*ILE*ITE*DYBU1*DYBU2.DYBL1*DYBL2. 

3 K*SMALLK»0MEG»XH.ND0UB»CPCPB*TITLE(8) * M8* del* alpha *ALPHAF 
COMMON  /COEFF/  A (99) .8 (99) *C (99) *D (99) .PHI ( 10000) 

COMMON  /GAMMA/  GAMTE1 *GAMTE*PGFF *GAMFF .GAMFFS 

PART=.5*( (X(ITE*1)-1.)*(GAMFF-GAMTE)/(X(IM1)-1.)*GAMTE) 

IJW=ITE*JM*JW 

PHI ( I JW ) =PHI ( I JW) -PART 

L=(ILE-2)*JM*JW 

PHIUB ( ILE-1 ) =PHI (L) 

PHIUB ( ITE* 1)=PHI (IJW)*2.*PART 
C COMPUTE  CP  LOWER  (B)  AND  CP  UPPER  (D) 

DO  297  I=IL£* ITE 
M= ( I - 1 ) * JM 
L=M* JW 

B(I)=-2.*(AX1 (I)*(PHI (L*JM)-PHI (L ) > ♦ AX2 ( I ) * (PHI (L ) -PHI <L- JM) ) ) 

D ( I ) =-2  • * ( AX  1 (I)* (PHIUB (1*1) -PHIUB ( I ) ) *AX2 ( I ) * (PHIUB ( I ) - 
1 PHIUB ( 1-1 ) ) ) 

297  CONTINUE 

PHI (IJW)=PHI (IJW) .PART 
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WRITE  (6.911)  NITERG 

WRITE  (6.903)  < D ( I > , I = ILE ♦ I TE > 

WRITE  (6.912)  NITERG 

WRITE  (6.903)  < B ( I ) . 1= ILE , I TE ) 

903  FORMAT  (10E13.5) 

911  FORMAT  (1H  ./.*  AT  I TERAT ION* 1 5*  SCALED  PRESSURE  COEFFICIENT.  UPPE 
1R  (ILE  TO  ITE)  =*) 

912  FORMAT  ( 1 H ./.*  AT  I TERAT ION* 15*  SCALED  PRESSURE  COEFFICIENT,  LOWE 
1R  (ILE  TO  ITE)  =*) 

RETURN 

END 

SUBROUTINE  TRI  (I) 

COMPLEX  PHIUB.B.D.PHI »P 
REAL  K.M8 

COMMON  /DELTA/  DX (99) .DY (99) ♦ AX1 (99) , AX2 (99) ,BX 1 (99) ,BX2 (99) , 

1 CX ( 99) , A Y 1 (99) .AY2(99) ,X(100) . Y ( 1 00 ) , FPU ( 99) .FPL (99) ,PHIUB(99) , 

2 I M , I M 1 ♦ JM ♦ JM 1 » JW. JWP1 ♦ JWM1 , ILE, ITE.DYBU1 .DYBU2.DYBL1 .DYBL2. 

3 K.SMALLK.OMEG.XH.NDOUB.CPCPB, TITLE (8) ♦ M8 ♦ DEL ♦ ALPHA , ALPHAF 
COMMON  /COEFF/  A (99) .B (99) ,C (99) ,D (99) .PHI ( 10000 ) 

DO  10  KK=3 ♦ JM 1 

J=JMl-KKO 

P=A(J-1)/B(J) 

B ( J-l ) =B ( J-l ) -P*C ( J) 

D ( J- 1 ) =D ( J-l ) -P*D ( J) 

10  CONTINUE 
M = < I — 1 ) * JM 
PHI (M.2)=D(2)/B(2) 

DO  20  J=3 » JM 1 
L=M*  J 

PHI (L)-(D(J) —PH I (L-1)*C(J) )/B(J) 

20  CONTINUE 
RETURN 
END 
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